1 1 ,< *>s
NfiVflL RESEARCH
LOGISTICS
QUflRTERLy
JUNE 1970
VOL. 17, NO. 2 fe
re
13 OCT 7 0%
OFFICE OF NAVAL RESEARCH
NAVSO P1278
NAVAL RESEARCH LOGISTICS QUARTERLY
EDITORS
Rear Admiral H. E. Eccles, USN (Retired)
The George Washington University
F. D. Rigby
Texas Technological College
S. M. Selig
Managing Editor
Office of Naval Research
Arlington, Va. 22217
O. Morgenstern
Princeton University
D. M. Gilford
U.S. Office of Education
ASSOCIATE EDITORS
R. Bellman, RAND Corporation
J. C. Busby, Jr., Captain, SC, USN (Retired)
W. W. Cooper, Carnegie Mellon University
J. G. Dean, Captain, SC, USN
G. Dyer, Vice Admiral, USN (Retired)
P. L. Folsom, Captain, USN (Retired)
M A. Geisler, RNAD Corporation
A. J. Hoffman, International Business
Machines Corporation
H. P. Jones, Commander, SC, USN (Retired)
S. Karlin, Stanford University
H. W. Kuhn, Princeton University
J. Laderman, Office of Naval Research
R. J. Lundegard, Office of Naval Research
W. H. Marlow, The George Washington University
B.J. McDonald, Office of Naval Research
R. E. McShane, Vice Admiral, USN (Retired)
W. F. Millson, Captain, SC, USN
H. D. Moore, Captain, SC, USN (Retired)
M. i. Rosenberg, Captain, USN (Retired)
D. Rosenblatt, National Bureau of Standards
J. V. Rosapepe, Commander, SC, USN (Retired)
T. L. Saaty, U.S. Arms Control and
Disarmament Agency
E. K. Scofield, Captain, SC, USN (Retired)
M. W. Shelly, University of Kansas
J. R. Simpson, Office of Naval Research
J. S. Skoczylas, Colonel, USMC
S. R. Smith, Naval Research Laboratory
H. Solomon, The George Washington University
I. Stakgold, Northwestern University
E. D. Stanley, Jr., Rear Admiral, USN (Retired)
C. Stein, Jr., Captain, SC, USN (Retired)
R. M. Thrall, University of Michigan
T. C. Varley, Office of Naval Research
C. B. Tompkins, University of California
J. F. Tynan, Commander, SC, USN (Retired)
J. D. Wilkes, Department of Defense, OASD (ISA)
The Naval Research Logistics Quarterly is devoted to the dissemination of scientific information in logistics and
will publish research and expository papers, including those in certain areas of mathematics, statistics, and economics,
relevant to the overall effort to improve the efficiency and effectiveness of logistics operations.
Information for Contributors is indicated on inside back cover.
The Naval Research Logistics Quarterly is published by the Office of Naval Research in the months of March, June,
September, and December and can be purchased from the Superintendent of Documents, U.S. Government Printing
Office, Washington, D.C. 20402. Subscription Price: S550 a year in the U.S. and Canada, S7.00 elsewhere. Cost of
individual issues may be obtained from the Superintendent of Documents.
The views and opinions expressed in this quarterly are those ofthe authors and not necessarily those of the Office
of Naval Research.
Issuance of this periodical approved in accordance with Department of the Navy Publications and Printing Regulations,
NAVEXOS P35
Permission has been granted to use the copyrighted material appearing in this publication.
MARKOV CHAIN ANALYSIS OF A SITUATION WHERE CANNIBALIZA
TION IS THE ONLY REPAIR ACTIVITY
Alan J. Rolfe
The Rand Corporation
Santa Monica, California
ABSTRACT
This paper considers a group of S identical aircraft, each of which is partitioned into
K parts which fail exponentially. The only way in which a failed aircraft can be repaired is by
cannibalizing its outofcommission parts from other failed aircraft. The evolution of the
number of good aircraft over time is governed by the transient behavior of an absorbing
Markov chain. We can therefore study this behavior by matrix multiplication although the
computational problem grows large for K 5= 3. Some numerical results and some approxi
mations are also provided.
1. INTRODUCTION
This paper presents a model of a system subject to cannibalization. The system is a simple special
case of the very general coherent structure considered by Hirsch, Meisner, and Boll (HMB) in [3].
It differs, however, from the model of HMB in its treatment of time.
In Sec. 2 the model is described in terms of a specific example. The reader can no doubt imagine a
number of other situations where the same model would be appropriate. We also make clear how our
model treats time differently from the model of HMB. Section 3 presents the Markov chain transition
matrix which characterizes our problem, while Sec. 4 explores the computational chore associated
with using the model. In Sec. 5 and 6, we exhibit several modest numerical examples and suggest
means of lessening the computational chore for larger problems. Section 7 contains a special technique
which is feasible for smallscale problems.
Although this paper does not completely and explicity solve the problem it poses, we hope that
the results given will provide quantitative insights into the effects of cannibalization which will lead
to further work in the field.
2. THE PROBLEM
We consider a group of S identical aircraft operating under the following conditions:
(a) Every aircraft includes K parts or assemblies. The timetofailure of the ith part is exponentially
distributed with distribution function, 1 — e~V. Futhermore, the K parts are assumed to be noncritical
items so that no matter how many of the K assemblies fail during a mission it is always possible for an
aircraft to return safely to base.
(b) A failed part can neither be repaired nor resupplied. However, there is sufficient repair capa
bility to cannibalize unfailed parts from other failed aircraft in an attempt to minimize the number of
outofcommission aircraft. We assume that initial supplies of spare parts are zero, although this
condition can be relaxed as we indicate later.
(c) The aircraft fly missions of fixed length T in which every "good" aircraft flies every mission.
We assume that there is enough time between consecutive mission takeoff times for all aircraft to
151
152 A. J. ROLFE
return from the previous mission and for all cannibalizations to be performed. We indicate below how
the "every good aircraft flies every mission" rule can be modified.
The basic quantity which the analysis below will provide is the expected number of good aircraft
after t missions, t=\, 2, . . ., denoted by E[A{t)]. This information will, for example, allow us to
measure the changes in force capability over time which cannibalization affords under the conditions
above.
It is important to distinguish the way in which our model differs in its treatment of time from that
of HMB. In our case, only those parts which happen to be in the air during a given period of time are
subject to risk. Good parts lying in partially cannibalized aircraft on the ground are not subject to
failure. In effect, each part carries its own clock which is turned on when the part is sent up into the
air and turned off when the part is on the ground. In the paper of HMB, however, every part that has
not previously failed is subject to risk. There is a single clock which is turned on at time 0, and while
it runs all parts are subject to failure. In Sec. 6 below we use the above fact to identify one of the
results of HMB as a lower bound for E[A(t)].
3. THE MARKOV TRANSITION MATRIX
If we let rii(t) be the number of copies of part i which survive t missions, then n(t) = (ni(t), . . .,
nA(O) represents the number of parts of each type which survive t missions. This is the basic state
vector of our process. As n.i{t) is an integer and ^ rii{t) ^ S for i=l, 2, . . ., K, the process has
(S + 1) K states for each t. However, for our purposes we can lump all states having one or more ni{t) =0
into a single "zero" state, thereby reducing the number of states to S K + 1. It is readily seen that the
process is Markovian and that its evolution from mission t to t+l can be described by a onestep
Markov chain transition matrix P of dimension S K + 1 by S K + 1. If the rows and columns are ordered
lexicographically, P will be a lower triangular matrix with 0's above the principal diagonal.
It is convenient to label row and column m by d m , where
d\ =(zero state);
<k =(1, . . ., 1);
d S K +1 =(S, . . .,S).
Then the matrix element py, where dt= (t'i, . . ., i/<), dj= (j\ . . .,Jk), and M,= min ii , is given by
i=\
K / M. \
(2a) f[ ' ,„ ) (e K i T ) j i i ' +M i(le K i T ) i iJi i
jJi \ji — ii + Mj/
if i,y>l and;, =S i, <;, + M„ for all /=1, . . .,K;
(2b) if i, j > 1 but ii < ji or i t > ji + Mi for at least one / = 1 , . . . , K;
(2c) lifi=j = l;
(2d) Oif £=1,j>1;
MARKOV CHAIN CANNIBALIZATION ANALYSIS 153
i
(2e) 1~2 Pidi£i>l,j=l
rf=2
These entries become apparent if we realize that Mi copies of each part are tested on each mission
and that each part is assumed to fail independently of other parts. Once the zero state is reached the
"experiment" ends as all aircraft are outofcommission.
4. THE COMPUTATIONAL CHORE
If we let ri(<) be the row vector of unconditional state probabilities after t missions, then from
Markov chain theory we know that
(3) n(*)=n(0)F, «=i,2, . . .,
where 11(0) = (0, . . .,0,1) as the process starts in state (S, S. . . . , S) and the states of II() are
ordered lexicographically as are those of P. Furthermore, P represents an absorbing Markov chain
[4] with d\ the absorbing state and all other states transient. From (3) we see that to go from Tl(t) to
Yl(t\\) we must multiply 11(0 , a 1 byS K + 1 matrix, byP, anS^f 1 byS A + 1 matrix. This operation re
quires roughly (S 2K /2) multiplications and the same number of additions. The major problem in using
the model is one of computational feasibility as (S 2A /2) can be a large number. For example, if S = 24
and K = 3, (S 2K /2)= 1.5 X 10 8 . Thus it appears that either S or K must be fairly severely restricted
if we are to perform the computation in an acceptable amount of time; the following sections describe
some numerical results which suggest ways of lessening the computational chore.
The first two moments of the number of good aircraft after t missions are given by
(4) E[A(t)]= S 2 nj(t)Mj
and j=1
(5) E[A(t)*] = S! 2 Uj(t)Ml
where
Mj= min {ji\dj={ji, . . ., >*)}•
1=1, . . .,K
Since Mj is independent of t for each;, the computation of E[A{t)] should proceed rapidly once
Yl(t) has been determined. As the sole absorbing state d\ of the Markov process has Mi = it is clear
from (4) that lim E[A(t)]=0; this is of course what one would expect given the norepair situation
being modeled. The computation of 11(f) would generally be iterated until E[A(t) ] was fairly close to
zero.
5. NUMERICAL RESULTS FOR K = 2
In order to simplify the computational task, we first limit ourselves to the case of two cannibalizable
parts per aircraft. There remain, then, two dimensions to the problem: S, the number of aircraft under
test, and g = {e~ x,T , e~ k2T ), the vector whose components are the probabilities of nonfailure of each
part on one mission of length T. Table 1 presents results for q fixed and S variable, while in Table 2
q is variable and S fixed. In both cases the quantity of interest is E[A(t)], the expected number of
good aircraft after t missions.
154
A. J. ROLFE
Table 1
\ s
t \.
Any S, no
cannibalization;
(0.72)'
2
5
10
13
1
0.720
0.734
0.762
0.777
0.783
2
0.518
0.557
0.600
0.622
0.628
3
0.373
0.431
0.480
0.500
0.504
4
0.268
0.336
0.384
0.402
0.405
5
0.193
0.263
0.308
0.323
0.326
6
0.139
0.206
0.247
0.260
0.261
7
0.100
0.161
0.198
0.209
0.210
8
0.072
0.125
0.159
0.168
0.169
9
0.054
0.097
0.128
0.135
0.136
10
0.037
0.075
0.102
0.108
0.109
11
0.027
0.058
0.082
0.087
0.087
12
0.019
0.044
0.066
0.070
0.070
13
0.014
0.034
0.053
0.056
0.056
14
0.010
0.026
0.042
0.045
0.045
15
0.007
0.019
0.033
0.036
0.036
(Entries are E[A(t)]IS; q = (0.9, 0.8). For the no cannibalization case
E[A(t)]/S = (0.72)' results from the fact that 0.72 is the probability that
one aircraft survives one mission and that we have a sequence of t bi
nomial trials.)
From Table 1 we see that for a fixed value of t, the ability to cannibalize produces an increase in
E[A(t)]IS which also increases with S. However, we note that as S increases the quantity E[A(t)]/S
Table 2
t \v
0.9, 0.8
0.8, 0.8
0.7, 0.7
0.6, 0.6
0.8, 0.5
0.5, 0.5
0.9, 0.4
0.5, 0.4
1
2
3
4
5
6
8
10
12
14
0.777
0.804
0.803
0.806
0.805
0.802
0.803
0.804
0.804
0.804
0.731
0.768
0.773
0.775
0.777
0.776
0.776
0.775
0.773
0.620
0.661
0.667
0.668
0.668
0.667
0.660
0.514
0.558
0.563
0.558
0.556
0.550
0.494
0.506
0.504
0.503
0.503
0.412
0.456
0.456
0.450
0.399
0.400
0.401
0.400
0.353
0.391
0.389
0.382
(Entries are E[A(t)]/E[A(t 1)]; S = 10.)
converges for values of S around 10. The entries for S — 13 are not much different from those forS = 10.
Hence it seems that we can learn as much about E[A(t)]/S by studying 13 aircraft as by studying 24 or
48 and the computational job is markedly easier. Of course, this conclusion is strictly valid only for
K — 2 and q = (0.9, 0.8). A question of interest (as yet unanswered) is, "At which value of S does con
vergence occur, given a particular K and q?"
In Table 2 we exhibit the behavior of successive values of E[A(t)] for a group of 10 aircraft and a
variety of gvectors. One fact is immediately apparent: For t 2= 2, E[A(t)]/E[A(t— 1)] =p, p a constant.
MARKOV CHAIN CANNIBALIZATION ANALYSIS 155
This is equivalent to saying £'[/i(f)]/S = p'" 1 , t Ss 2, which implies
(6) E[A(t)]IS = p ( P y\t^l,
where
pomEtA(l)]IS.
This provides an immediate analogy to the nocannibalization case, where we know that E[A{t)]/S = A.',
«Sl, with A = <7i<7 2 where q = (qi , q 2 ) . (It was found that for S = 5 the geometric behavior was less stable
than that characterized by Table 2 for S= 10).
The next section will suggest several possible ways around the computational problem.
6. APPROXIMATIONS FOR K>2
In Sec. 5 we indicated that for K= 2, any S 3= 10 and any q we could approximate the behavior of
E[A(t)]/S by studying that of the case S— 10. A typical computer run on the 360/65 for S= 10 and
q= (0.9, 0.8) took approximately 7 sec. If K were increased to 3, the information in Sec. 4 would
suggest a run time of roughly 10 min., still a feasible amount of time. Hence for K= 2,3 we can study
a group of 10 aircraft and be fairly confident that the E[A (t)]IS so obtained closely approximates that
for any larger value of S and the same q.
For K 25 4, however, studying even 10 aircraft becomes prohibitively costly as each time K is
increased by 1, computation times are increased by a factor of 100. Thus, although the geometric
behavior exhibited by E[A(t)]/S in (6) would probably still hold, the problem of determining p and p
could no longer be handled by numerical means.
An approximation for E[A(t)] is given by (7), with q= {q u ■ ■ • , <7a ) and a=mjn </,,
(7, 5f*^,„i
Relation (7) is based on the imprecise reasoning that A(t) is equal to the number of survivors of the
parttype for which qi is minimized.
The numbers in Table 2 suggest that the RHS of (7) is often an upper bound for E[A(t)]IS rather
than an approximation; however, column 5 of that table shows that this is not always so. For t=l, the
relationship E[A( 1)]/S =S a follows readily from min Xj^Xi, i=l, . . . , K,
i = 1 K
K
A lower bound for E[A(t)]/S is given by (8), with 8 = r[ qi,
The referee suggested the following better lower bound, which can be derived from formulae 4.3
and 4.9 of [3]. This result can also be derived directly without reference to [3] and is based on the
supposition that every surviving part be tested on every mission.
(9) E[A(t)]
The referee provided the results shown in Table 3 whose entries are the lower bound, (9), for the case
shown in Table 1. Expression (9) appears to be an excellent lower bound in this instance.
156
A. J. ROLFE
Table 3
\ s
t \^
2
5
10
13
1
0.734
0.759
0.777
0.782
2
0.554
0.598
0.620
0.626
3
0.423
0.473 ■
0.495
0.501
4
0.323
0.374
0.396
0.401
5
0.247
0.296
0.316
0.320
6
0.187
0.233
0.252
0.256
7
0.142
0.184
0.201
0.204
8
0.106
0.144
0.160
0.163
9
0.080
0.113
0.127
0.130
10
0.059
0.088
0.101
0.103
11
0.044
0.068
0.080
0.082
12
0.032
0.052
0.063
0.065
13
0.024
0.040
0.049
0.052
14
0.017
0.030
0.039
0.041
15
0.013
0.023
0.030
0.032
7. A RECURSIVE RELATION FOR THE GENERATING FUNCTION OF THE STATE
PRORABILITIES
The fact that the Pmatrix, as given by expression (2), is lower triangular allows us to develop a
recursive relation which, at least in theory, produces a closed form expression for the generating func
tion of the state probability vector, Yl{t). Inverting these generating functions presents a barrier if
S and Kare large; however, we present the method as it might sometimes be desirable to study a problem
with small S and K.
Setting N=S K +l, we let {ptj, i, j=l, . . ., N} be the elements of the Pmatrix. If II(t) = (IIi(f),
. . ., Il/v(t)) is the vector of state probabilities, we define the generating functions
(10) Ui(Z) = £ Yli(t)Z', i = l, 2, . . ., N.
t =
The basic Markov chain relation is
(11)
nu) = n{t\)P.
Using (10), (11), and the fact that pjj = 0,j> i, we can derive
U N (Z) =p NN ZU N (Z) + 1
This is equivalent to
(12a)
(12b)
Yl i (Z) = ^ Pji Zn j (Z), i=Nl,N2, . . ., 1.
n N {Z) = \j{\p NN z)
n,(Z) = (z 2 p ji n j (Z))i(ip ji z), i=Ni,N2, . . ., 1.
j = i + l
MARKOV CHAIN CANNIBALIZATION ANALYSIS
157
From (12a) we know at once that
Us(t) = (p NN y, t^O.
Then (12b) can be solved recursively as IT(Z) depends only upon flj(Z), j>i. The difficulty in a
large problem is that inverting Yli(Z) is not trivial as it was in (12a) because the generating function will
have a large number of multiple roots. There do exist both exact and approximate methods for such
problems, some of which are found in [2, Chap. XI]. We present a short numerical example to illustrate
the technique:
S =2
K = 2
q= (0.5, 0.4)
From expression (2), we find that
P =
1
0.8
0.2
0.5
0.3
0.2
0.6
0.2
0.2
_0.52
0.24
0.08
0.12
0.04
Then expression (12) produces
n 5 (Z) = i/uo.04Z)
fl 4 (Z) = (0.12Z)/(1 0.2Z) (1 0.04Z)
n 3 (Z) = (0.8Z)/(1 0.2Z) (1 0.04Z)
n 2 (Z) = (0.24Z)/(l0.2Z) 2 (l0.04Z)
n I (Z) = (0.0016ZM5 + Z)(65Z)/(lZ)(l0.2Z) 2 U0.04Z).
If we use methods in [1] to invert these generating functions and then use the relation
we find
E[AU)] = o{n l u)} + i{n 2 (t) + iu(t) + n 4 (t)} + 2{Ti i (t)}.
^(*)]=«(.2V+(.2)'+(.04)',«s»0.
Thus we have found a closed form expression for E[A(t)].
8. VARIATIONS OF THE MODEL
Several of the assumptions made in Sec. 2 can be relaxed without materially affecting the develop
ment outlined above. The only change is in the Markov transition matrix P:
158 A. J. ROLFE
(a) There may be an initial supply, s,, of spares for part i, i= 1, 2, . . . , K. In this case the possible
states of part i are 0, 1, . . .,5, + S and the Pmatrix has
f[( Si + S) + l
rows and columns. This would make the computational task mentioned in Sec. 4 more difficult;
(b) The assumption that "every good aircraft flies" every mission" can be replaced by the re
quirement that only Q aircraft, if available, fly every mission. This would be accomplished by sub
stituting Q for Mi in expression (2a) whenever Q < Mi. If this were done, the Pmatrix would be made
more sparse while E[A(t)] would approach zero less rapidly because the Markov process would take
more time to work its way from the initial state to the absorbing state.
9. ACKNOWLEDGMENTS
The author would like to thank Dr. D. M. Landi of The Rand Corporation for his part in bringing
the model to the author's attention. The author is especially indebted to an anonymous referee who
supplied many useful comments, all of which are incorporated in this paper.
REFERENCES
[1] D'Azzo,J. J. and C. H. Houpis, Feedback Control System Analysis and Synthesis (McGrawHill, New York, 1960).
[2] Feller, W., An Introduction to Probability Theory and Its Applications (Wiley, New York, 1957), Vol. 1.
[3] Hirsch, W., M. Meisner, and C. Boll, "Cannibalization in Multicomponent Systems and the Theory of Reliability," Na'v.
Res. Log. Quart. 7.5, 335359 (1968).
[4] Kemeny, J. G. and J. L. Snell, Finite Markov Chains (Van Nostrand, Princeton, New Jersey, 1960).
SOME ESTIMATES OF RELIABILITY USING INTERFERENCE
THEORY
M. Mazumdar
Westinghouse Ktscarch Laboratories .
Pittsburgh, Pennsylvania 15235
ABSTRACT
According to "interference theory" of reliability, a component fails if the maximum
stress exceeds the component's strength. Assuming that both these quantities are random
and their distributions are normal, we obtain in this paper some point and interval estimates
of reliability when the stress distribution is known and a few observations exist on com
ponent strengths.
1. INTRODUCTION
In recent years, many writers have attempted to analyze the reliability of a component by the
application of probability arguments to a certain physical model of failure. According to this model
(which has come to be recently known in the literature as "interference theory"), a component fails
if at any moment the applied stress (or load) is greater than its strength or resistance. The stress is a
function of the environment to which the component is subjected, and its value at any point of time
is considered to be a random variable. The strength of the component is measured by the stress re
quired to cause failure. Strength depends on material properties, manufacturing procedures, and so on.
If the components under question are mass produced and their selection in a given system is assumed
to be made at random, then the strength should also be considefed a random variable. The reliability
of a component during a given period [0, T\ is taken to be the probability that its strength exceeds the
maximum stress during the entire interval. The above formulation of the reliability problem was
first proposed by Birnbaum [1], who also gave a method for estimating the related probability. In
recent years, the model has found an increasing number of applications in many different areas,
especially in the structural and aircraft industries.
The stress to which a component is subject in its operating environment can often be predicted
using the available technological knowledge about the mechanical, thermal and other relevant condi
tions of the system, and the way they interact with each other. In many cases mathematical models
are available which allow one to compute deterministically the values of "stress" at different points of
time corresponding to a given set of input values describing the initial conditions of the system. For
example, as Church and Harris [2] have pointed out, in a missile flight these initial values may cor
respond to the propulsive force, angles of elevation, changes in atmospheric conditions, etc. Since
these input values have statistical distributions associated with them, the calculated "stress" values
will also be randomly distributed. Moreover, since it will be possible to generate a practically unlimited
number of such stress values using simulation methods by varying the input values in accordance with
their respective distributions and hence to reconstruct the stress distribution as accurately as desired,
we might assume that the distribution of the maximum stress is known. The strength distribution can
not, however, usually be computed from such a priori considerations. It can only be estimated using
159
160 M. MAZUMDAR
statistical methods from the results of tests carried out for this purpose. The number of such tests is
likely to be small in most cases due to practical limitations.
Under, the assumption that the distributions for strength and the maximum stress are both normal,
we derive in this paper some point and interval estimates of reliability using "interference theory,"
when as above the stress distribution is known and the strength distribution is unknown, but a few
observations on the latter are available. The point estimates are given in section 2, and the interval
estimates are derived in section 3.
Owen, Craswell, and Hanson [7] considered a similar problem and gave confidence intervals for
the reliability when both stress and strength are (a) jointly bivariate normally distributed and observa
tions are in pairs, or (b) when they are independently normally distributed, but they have a common
(but unknown) variance. Most recently, Church and Harris [2] considered the same problem and
provided a procedure for obtaining confidence intervals, which they compare with that of Govindarajalu
[3]. Their procedure, although empirically demonstrated to be superior to that of [3], is, however,
inexact since it uses the asymptotic normal approximation of a given statistic and requires the substi
tution of the population mean and standard deviations by their observed sample values. Our procedure
for obtaining the confidence interval is based on that given in [7] and, though exact in a certain fashion,
is admittedly inefficient. Some comparisons of this procedure with that of Church and Harris [2] are
given in section 3. Sampling experiments show that the average length of the interval associated with
the Church and Harris procedure is smaller than that of our procedure although it appears that, in
general, it fails to cover the true parameter value more often.
2. POINT ESTIMATES
Let X denote the maximum stress to which the component is subject over a given period of time.
We assume that X is normally distributed with mean /jli and variance cr'f, and that both these values
are known. Let Y denote its strength. We assume that Y is also normally distributed with unknown
mean /xz and variance erf. In what follows, we shall first assume that cr is known and then that it is
unknown. No doubt the latter is the more realistic case. In terms of our definitions we readily see the
reliability to be
(1) p = Pr{YX>0}
= <b( ^^ )
\ Vof + of/ '
where
®(t)=—^= (' e x *' 2 dx.
Let Yi,Yz, . . . , Y m denote m mutually independent observations on the strength of the components.
These values are obtained from bench tests. On the basis of these observations, we obtain below the
minimum variance unbiased (MVU) and the maximum likelihood estimates of p.
2.1 Case 1: cr known
Let Z be a normal random variable with mean and variance erf, and let
(2) / =
1
if
Y x + Z> /x,
if
Yi + Z^m
RELIABILITY BY USING INTERFERENCE THEORY 161
Then
(3) E{I}=Pr{Y l +Z> f x 1 }
=P
It therefore follows from BlackwellRao and LehmannScheffe theorems (See Chap. 3. [5]) that
_ _ »'
the conditional expectation E{I\Y} is the MVUestimate of p, where F=VF;/m is the sufficient
i= 1
statistic. Since moreover the conditional distribution of Y t + Z given F is normal with
(4) £{F,+ZF} = F,
Var{F 1 +ZF} = cr2 + o('l) ,
the MVUestimate of p in this case is given by
(5) p, = E{h\Y)
= Pr{Y i + Z>p. 1 \Y}
= *( ^ ).
VVof + ofdl/m)/
The maximum likelihood estimate of p, on the other hand, is given by
It is seen that as m becomes large, the two estimates become equivalent.
2.2 Case 2: erf unknown
In this case the sufficient statistics is the pair (F, S.J) where
»' _
Si=2 (F,Y)*.
Therefore, using arguments identical to those of Case 1, it follows that in this case, E{I/Y, S'i} is the
MVUestimate of p. From [6] we find the conditional density of Y\ given F and Sf to be
' m \ 2 ) 1 f (YyV m ] 2
m — , „/m — 2\ S» 1 \ S> J m —
otherwise.
162 M. MAZUMDAR
Therefore the MVUestimate of p in this case is
(8) p n = E{I\Y,Si}
= Pr{ Y x +Z >/xi  Y,S%}
V m
I i Y—y / m
By use of the transformation 2 =  +  —  — \  the integral on the right hand side of (8) reduces to
z z 02 v m — I
Jo r
(9) />// = I ' . r( o m , 2) . zf"'"/ 2 (lz) (, " 4)/2 ^(« + j82)ffe,
'm — Z\ v (m — Z x
where a =(s.,J— 1 Y—pn )/cr,,
m
and /3 :
2S 2 \m —
<j\ \ m
It does not appear that this integral can be evaluated in a closed form; however, it can be evalu
ated very easily for particular examples by the use of numerical integration.
Thf* maximum likelihood estimate forp in this case is
(10) P// = <*>'
where
Vaf+4
S?
ss =
2 mr
pu is biased, but might be quite efficient when m is reasonably large. Also, it can be shown that for
large values of m, the MVUestimate pu and the maximum likelihood estimate/?// are equivalent.
3. INTERVAL ESTIMATES
As in the discussions in the preceding sections, we give the interval estimates of reliability for two
separate cases (a) when cr\ is known, and (b) when o" 2 is unknown. We derive here only the lower
confidence limit because of its importance in the present application. An extension of our procedure
to twosided intervals will be obvious.
3.1 Case I: o~ 2 known
Let K a be such that <$>(K a ) = a. Since Y is normal with mean p., and variance cr 2 /m, it follows that
,11) Pr{ V "<r*> «*.}!„,
f fii — fiz Ui — Y—K a <T 2 lVm} .
or rr \ , =£ \ = 1 — a,
I V o? + o 2 V rr? + 02 J
'oi + 0 2
RELIABILITY BY USING INTERFERENCE THEORY 163
Pr
k <^1 h
Thus the 100 (1 — a) percent lower confidence limit for p in Case I is
(12) p,,dla)=<t>l ./—; )•
\ Vof + crj '
3.2 Case II: &\ unknown
The procedure suggested by Church and Harris [2] for obtaining the interval estimate in this case
is based on the statistic
v ? ^
^ai + sl
They have shown that V is asymptotically normal with mean (ju, 2 — Pi )/ V^f^ and variance
(/Lt 2 — (JLi) 2 <ri
OV —
m 2(ml)(o'f + o2) 2
1 ■ "2>
Using normal approximation for the distribution of V and upon replacing the unknown values of fx 2
and <t\ by their respective sample estimates Y ands*, they obtain (the symbol = denotes approximately
equal to)
(13) Pr { / *IZ^ > V+ K a CTv] = l«,
I Vo2 + O"! J
or Pr{p>4>(K+^„ov)} = 1a,
or Pr{p><D(r+^)}=la,
where
&l ai + ,\
■+
m 2(m\)((T 2 + siy
Since this procedure involves two approximations both of which may not be close in the nonasymptotic
case, the above results may not be satisfactory in small samples.
We next suggest a procedure based on results given in section 4 of [7]. This procedure is "exact,"
but is admittedly inefficient because it does not make use of all the available information; however, the
use of sampling adopted in this method seems to be the only way we can eliminate the nuisance param
eter &\ without introducing any approximations, and we believe that the properties of the procedure
might therefore be of interest. It seems possible that ad hoc procedures can be developed which use
this sampling concept, but are more efficient than the procedure described here. Results obtained from
simulation experiments given below seem to indicate that our procedure covers, in general, the true
parameter value more often than the Church and Harris procedure, although the average interval
obtained from the latter procedure is significantly shorter.
According to our procedure, we first draw, using random number tables, a sample of m independent
observations form a normal distribution with mean and variance cr\. Having taken these m observations
164 M. MAZUMDAR
we consider the transformed observations
It is clear that
(14)
and
Now let
(15)
and
Then
(16)
Wi = Yi + Zi, i = l, 2, . . ., to.
E{Wi}=p 2
Vai{Wi}=(r\ + trl
w=^ —
TO
£ {Wiwy
m^\
Sw
Sw
Vof + of Vo, 2 + o 2 2 J/ VqHd/
is distributed as a noncentral < with degrees of freedom = to— 1 and the noncentrality parameter
8— Vto(ju,i — /A2)/ Vof* + o" . Let the left hand side of (16) be denoted by T,,, 1. It then follows from the
formulas given in [4] that
(17) Pr{ ^ (/"**»> >S(ml < r ffl  t , q) j = lq,
( Vo 2 + crl J
where from Eq. (33) of [4], 8(to— 1, T m u ot) is defined as
(18) 1 8(to — 1, r m _i, a) =T m i — \(m— 1, r m _i, a)
1 +
71 ,
2(to1)
1/2
and k(m— 1, r m _i, a) can be obtained from Table IV of Johnson and Welch [4]. From (17) and (18)
it now follows that
5(to — 1, T m u a)
(19)
Pr
{p>^(
TO
= la,
so that a lower (1 — a) percent limit forp following this procedure is given by <J>(S(to— 1, T m  U a)/Vm).
RELIABILITY BY USING INTERFERENCE THEORY
165
We now exhibit the results obtained from 200 artificial sampling experiments comparing the
above two procedures using the same set of random numbers. In Table 1 the Church and Harris
procedure is denoted as Procedure A and the method suggested in this paper as Procedure B. The
particular sample sizes shown in the table were chosen because of the comparative ease with which
the values of 8(m — 1, T m i, a) can be evaluated from the existing tables.
Table 1. Comparison of Confidence Interval Procedures
0*i=
0, o, = l; 1
a = 0.95)
p
pz
0"2
m
Number of onesided con
fidence intervals out of 200
which covered true param
eter value
Average length of
200 confidence
intervals
Proc. A
Proc. B
Proc. A
Proc. B
0.90
1.43287
0.5
10
184
191
0.15
0.28
37
183
195
0.12
0.19
1.81246
1.0
10
171
187
0.19
0.26
37
185
193
0.14
0.19
4.05278
3.0
10
187
183
0.25
0.27
37
187
191
0.17
0.18
0.95
1.83905
0.5
10
181
188
0.087
0.19
37
188
190
0.067
0.10
2.32624
1.0
10
184
194
0.12
0.21
37
184
190
0.08
0.11
5.20163
3.0
10
187
189
0.18
0.20
37
185
189
0.10
0.11
0.99
2.60055
0.5
10
178
187
0.022
0.095
37
186
193
0.015
0.038
3.28946
1.0
10
182
194
0.045
0.097
37
183
189
0.022
0.039
7.35546
3.0
10
190
193
0.093
0.11
37
192
188
0.034
0.039
REFERENCES
[1] Birnbaum, Z. W., "On a Use of the MannWhitney Statistic," Proceedings of the Third Berkeley Symposium on Mathematical
Statistics and Probability (Uni\ersity of California Press, Berkeley, Calif., 1956), Vol. 1, pp. 1317.
[2] Church, J. D. and Bernard Harris, "The Estimation of Reliability from StressStrength Relationships," Technometrics, 12,
4954 (1970).
[3] Govindarajalu, Z., "DistributionFree Confidence Bounds for P{X < Y}," Annals of the Institute of Statistical Mathematics,
20, 229238 (1968).
[4] Johnson, N. L. and B. L. Welch, "Applications of the NonCentral t Distribution," Biometrika, Vol. .31, 362389 (1940).
[5] Lehmann, E., Notes on the Theory of Estimation, Mimeo Notes (University of California Press, Berkeley, Calif.) (Reprinted
1962).
[6] Lieberman, C. J. and G. J. Resnikoff, "Sampling Plans for Inspection by Variables," Journal of the American Statistical
Association, 50, 457516 (1955).
[7] Owen, D. B., K. J. Craswell, and D. L. Hanson, "NonParametric Upper Confidence Bounds for P{Y < X) and Confidence
Limits for P{Y <X} When X and Y Are Normal," Journal of the American Statistical Association, 59, 906924 (1964).
THE OPTIMAL BURNIN TESTING OF REPAIRABLE EQUIPMENT*
John M. Cozzolino
Department of Statistics and Operations Research
University of Pennsylvania
Philadelphia, Pennsylvania
ABSTRACT
The subject of this paper is the utilization of the "infant mortality" or decreasing failure
rate effect to improve the reliability of repairable devices. Decreasing failure rate implies
the possibility that devices which exhibit it can be improved by "burnin testing" of each unit.
Such a test serves to accumulate operating time while shielded from the full costs and con
sequences of failure. A general formulation of the burnin test decision for repairable devices
is presented and some special cases are solved. A class of models, indexed by the degree of
partial replacement present in the repair process, is considered and numerical results for
the optimal policy are given for several members of that class. A comparison of those results
reveals the profitability of testing increases with the complexity of the repairable device.
1. INTRODUCTION
The subject of this paper is the utilization of the "infant mortality" effect to improve the reliability
of repairable devices. The term infant mortality denotes a decrease with age of the conditional proba
bility of failure of a device and, consequently, an improvement in reliability with time as the device
operates. This phenomenon is widely observed. For example, decreasing failure rate has been found
to occur in electronic systems and components [12], biological systems [5], social systems such as
business enterprises [20], and also, in mechanical systems. In the case of repairable equipments de
creasing failure rate is often observed after the occurrence of failure and repair.
Decreasing failure rate implies the possibility that devices which exhibit it can be improved by
"burnin testing" of each unit. The term burnin testing is taken to mean test stand operation under
environmental conditions as close as possible to true operating conditions. The basic purpose of this
test is to accumulate operating time while shielded from the full consequences and costs of failure.
The measurement of reliability, while not a primary goal, may be a partial byproduct.
The main theme of this paper is a general formulation of the burnin test decision problem for
repairable equipment. A continuous stage dynamic program is used to treat this sequential decision
problem and some solutions are given for special cases. Space does not permit a full discussion of the
class of models for which this decision structure applies; however, one broad family of models, indexed
by the degree of partial replacement present in the repair process, is presented. Numerical results
for the optimal test policy are given for several members of that family, and a comparison of those
results reveals an increase in the profitability of testing as a function of the complexity of the repairable
device.
*The work herein reported was supported by the Office of Naval Research under Contract Nonr 1841 (87); NR 042230 while
the author was at the Operations Research Center at the Massachusetts Institute of Technology. The computation reported here
was done in part at the Computation Center at Massachusetts Institute of Technology.
167
168 J. M. COZZOLINO
2. LITERATURE REVIEW
The distribution theoretic conditions under which the expected future life of a device increases
with its age are shown by Watson and Wells [27]. Based upon the work of Knight [18], one stated
result is that the expected future life of an aged device always exceeds that of a new one if the device
has a monotonically decreasing failure rate.
Several authors have investigated the burnin test problem for the case of unrepairable equipments.
In that problem a failed device is considered worthless. The quantities investigated are the burnin
time required to achieve a specified expected future life of the surviving units. Watson and Wells [27]
discuss these quantities for several failure distributions. A nonparametric method is considered in
a paper by Lawrence [19]. Assuming only a monotonically decreasing failure rate, he derives bounds
upon the burnin time required to achieve a specified expected future life.
Recently, interest has developed in the burnin testing of repairable devices. This potentially
important problem requires some explicit assumptions about the nature of the relevant failure and
repair processes. This necessity follows from the need to assess the relative worth of the repaired
device. Kooharian and Young [28] consider a useful class of reliability models for repairable complex
devices. The outstanding feature which they prove for this class is a special simplified form of the
burnin test decision policy. A similar model for the "debugging period" is given by Bazovsky [4].
Concern for the long run reliability improvement of complex systems is evidenced by a paper by
Barlow and Scheuer [3] entitled "Reliability Growth During a Development and Testing Program."
They attributed the improvement mainly to engineering changes.
In a report upon which this paper is based [7], the author presents some models of decreasing
failure rate processes incorporating explicit consideration of the repair process. This union of failure
and repair models appears to be a fruitful basis for the discussion of burnin testing in view of the
interaction between the two processes. A complete picture of the life of the device is given in terms
of a stochastic process, usually with nonindependent interfailure times.
3. ELEMENTS OF THE PROBLEM
The Unit Test
A unit test starts with the testing of a new unit and ends with the output of one tested unit. If a
failure occurs during test, the test is interrupted, the unit is repaired, and the test then resumes. The
time elapsed during repair is not counted.
The Failure Process for New Units
Given a population of indistinguishable new devices, a specified and constant environment in
which they are intended to operate, and a set of conditions which define failure, a distribution function
Fo(tl6) with parameters 6 and density fo(t/6) is utilized t6 represent the uncertainty which exists re
garding t, the duration of operating time prior to the occurrence of failure for any one of the devices.
The failure rate of that density, also known as the hazard rate, conditional probability of instantaneous
failure, or force of mortality, is a very important function for determining the reliability characteristics
of the device. For a distribution satisfying Fo(0~/6) —0, the failure rate, a function of time, is defined as
It will be assumed throughout thatfo(t/6) is a continuous function of t and that Fo(O~/0) = 0.
OPTIMAL BURNIN TESTING « 169
Density of Future Time
The original density, /o, represents the process only before testing begins. Suppose that at some
time, r, after the test starts no failure has been observed in (0, t). Let t = future time to failure. Then
the total time to failure is (£+ r). Of this total time to failure the only random part remaining is t. The
relevant density function is that of future time to failure. This renormalized form of the density will be
denoted here by
f mr) M,+Tle) 
[lF„(r/«)]
While this density function usually represents a larger family of densities than does/oU/0), its failure
rate is the same except displaced to the left by t; no(t + Tl8).
The Repair Process
We consider that a failure during test will be followed by a repair of the failed unit. The rule which
assigns a new density of time to next failure for the repaired unit will be called the "repair rule."
In many analysis it is assumed that the repaired unit has the same density of time to next failure
as the new units. This practice ignores the information available from the failure and from the past
history of the unit. Since it is possible to distinguish between new and repaired units we realize that
it may be valuable to do so. This is particularly relevant for burnin testing models since the purpose
of testing is to create a difference between new and aged units.
The only restriction which we shall place upon the new density is that it be from the same family
of density functions as is f{t/8, t). Hence we allow t to have the meaning of a parameter of the density
rather than strictly the time since last failure.
The repair rule then is a transformation which determines which density of that family applies.
The most general transformation to be considered here is one that allows the new parameter values
to be a function of the parameter values which applied just prior to the failure. Hence the new parameter
values are
d R =o R (o,T)
and
tr = t h {6, t),
where (8, r) are the parameter values just prior to failure.
These functions will be required to have the property that, when applied to any point in the space
of allowed parameter values, they will produce a point which is also inside that space. Also, their con
tinuity in t will be required.
The repair rule results from an explicit model of the repair process. Examples of this are given
in [71
The Reliability Value Function
The reliability value function specifies the value that the tested unit would have if testing were
terminated. It gives the value as a function of any possible test position (6, t). It will be denoted by:
U(8, t) = the values of a unit when testing is terminated and a time when the state of process is {8, t).
This class of functions, U(8, r), includes all expectations of the future behavior of the unit, since
it includes all the parameters of the density function of future life.
170 J M. COZZOLINO
The Cost of Testing
The last remaining element of the problem is the cost of testing. In this we seek to represent
the actual cost of using the test facilities. The cost of testing will be considered to consist of the fol
lowing two components:
Ci = The cost of repair of a unit which fails during test. When several independent failure modes are
considered this definition can be expanded to include a different cost for each failure type.
C2 — The cost per unit test time.
4. OPTIMIZATION BY DYNAMIC PROGRAMMING
The profit from a unit test is the terminal value (given by the reliability value function) minus
the test costs incurred during the testing process. At an intermediate time during the test the optimal
test policy does not depend upon the costs already incurred since they are fixed or "sunk" costs.
The optimal policy is that which maximizes the expected profit of the remaining portion of the unit test.
The decisionmaking process is sequential; the decision to stop or continue is remade after each
observation of the partial test results. Since partial test results are revealed continuously in time each
time point is, in fact, a decision point. The basic repeated unit of the decision tree is shown in Figure 1.
/ /
, DEC IS ION NODE
STOP TEST FAILURE
NODE
•O CHANCE NODE
NO FAILURE
DECISION NODE
FIGURE 1. The basic unit of the Decision Tree
We define the optimal value function as
V(6, t) =the expected value of the remaining portion of the unit test when the current state is (6, t)
and when all future test decisions are optimal.
Let us consider a small positive time interval, h, starting at a time when the state of the process is
(6, t). The optimal value function, V(6, t), can be expressed in terms of the optimal decision at the
start of the interval, the events which can occur during it, and the optimal value function at the possible
resulting states. Let us consider the possibility of one event occurring in (t, t + A) at time r + e, where
e < h.
Since there are two possible decisions the optimal' value function is equal to the maximum of
two terms. The equation is
V(0, T) = Max{tf(0, t),C 2 A +
F (T + hlO)Fo(Tldy
1Fo(t/0)
[d + VlOniO, r+e), t h (6, r+e))]
lF (T+h/0)
+
lFo(rld)
V(d,T+h)
The term V(d, t) can be subtracted from both terms in the Max bracket. The resulting expression can
be written as
0=Max (J(0, t)V(0, t), h \C 2
,t),a{<
OPTIMAL BURNIN TESTING 171
V(d, T+h)~V(d,T)
+
F (T+hld)~Fo(Tld)
h(lF o (Tl0))
[C l + V(6 R {0,T + e),T R {d,T + e))V(0,T+h)]\\.
Let us now observe that we have the form
O=Max{£/(0, t)V(6, t), h£(h)},
where £(h) represents the interior of the second term. It should be noticed that the value of A has
no effect through the term £(h) since when hi;{h) <0 the function V(6, t) is determined by the
equation 0= U(0, t) — V(6, t). In fact, for any positive h,
h£(h)<0 <^f(A)<0,Af(A)=0<^f(A)=0.
Therefore, we can remove h and write the equation as
O = Max{£/(0, t)V{6, t), £(h)}.
The interval, h, can be considered as an inspection interval with the understanding that if a failure
occurs anytime during the interval, its occurrence will not be detected until time, t+ h, at which time
the operation of the device will be resumed. If detection is not an important aspect of the problem,
however, the limiting form of the equation is more useful. Provided only that the functions /o(t/0),
6h(6, t), and tr{6, t) are continuous functions of t everywhere in their range of definition, we can
now let h approach 0. In the limit, we obtain the equation
0= Maxjf/(0, t)  V(6, r) , dV{ ^ T)  «o(t/0) [V(6, t)  V(d«, r«) + C,] cX
It is helpful to visualize the testing process as a point moving with time in the (0, t) space. As
testing proceeds, r (effective age of the unit) increases and remains constant at nonfailure points
in time. When a failure occurs, both and r normally jump to a new value given by (0«, Tr). The
solution of the dynamic program divides the (0, t) space into two types of regions:
1. Continuation region where
(CI) U(6, t)V(6, t) «0, and
(C2) ^ ( ^ t) /xo(t/0)[F(0, T)r(0 ft ,T H ) + C,3C 2 = O.
2. Termination region where
(21) U(d, r)V(d, t)=0, and
(T2) d ^ T)  M o(T/0) [^(0, r)  F(0«, r«) + C,]  C 2 ^ 0.
These regions are separated by a boundary, t*(0), which represents the decision to terminate
testing. A simple illustration of a process trajectory in state space is given in Figure 2.
172
J. M. COZZOLINO
STATE SPACE (ONE 9 PARAMETER)
•(e)
REPAIR JUMP
TERMINATION
BOUNDARY
o Termination Point
X Failure Occupied
FIGURE 2. A process trajectory in state space
The following theorem helps to characterize the boundary in terms of the properties of U(6, r).
THEOREM: For each 6, the optimal policy, t*(0), is restricted to the following possibilities:
a) t* = or oo.
b) t* = t such that the derivative, 2 , is discontinuous at t.
St
, * . , dU(6, t) .
c) t =t such that 2 is continuous at t,
dr
d U(0, t) 8V(6,t) , a 2 t/(6>, t) d 2 V{d, t)
, and — s= — at r.
8t
dr
dr 2
3t 2
PROOF: Assume that, for given 6, V{6, t) and t*(^) are the optimal value function and decision
boundary, and that
a) 0<t*(6>) <oo
b) U(6, t) has continuous first derivative,  1 at t*.
OT
At a point {6, r* — e) (for e > 0) in the neighborhood to the left of t*, both V(d, t) and U{6, r) and
their derivatives exist by assumption. Condition CI gives
U(6,T*e)V(d,T*e) ^0.
Expansion about the point t = t* gives
U(6, T*e)  V(6, T*e)=U(d, r*)  V(d, r*) e
dU(d,r*ye) dV(6,T*ye)
dr St
for some < y < 1. Conditions CI and T\ then imply the following condition in terms of the lefthand
derivatives
dV{e,T*) ^ dU(e,T*)
dr dr
The lefthand derivative also satisfies Eq. C2
ay (e, t*)
dr
fjLo(T*/e)[V(e, T *)  v(o R , t«) +c] c 2 =o.
In a neighborhood to the right of r* the conditions 71 and T2 give
dU+(d, T*)
dl
■ti (T*id) [U(e, T )  v(o K , th) +C] c 2 ^ 0.
OPTIMAL BURNIN TESTING
Subtracting the latter two conditions and using 71, we have
dV(e,T) dU+(e, t)
173
dr
dT
Assumption C ensures that both onesided derivatives of U are equal; hence we have the existence
of Ps derivative at t* and the equality
dV(e, r) dU(e, t)
dT
dr
atr = r.
Furthermore, if the functions U and V have second derivatives at t* then a Taylor series expansion
using both derivatives and the above equation yields
d 2 V(d, r) d 2 U(d, t)
dT 2
dT 2
at r = r*.
Hence part Q of the theorem is proved. The possibilities excluded in assumptions a and b are parts
a and b of the theorem.
The differential equation is a necessary condition upon the optimal value function. Its form and
complexity will depend upon the nature of the 6 space and repair functions. In the simplest cases the
complete family of solutions of the differential equation can be obtained analytically. The parameters
V(8,r) for fixed 8,

1
1 ___
L_ " V(8,t) = U(8,t)
V(8,r> _______ ;  ; ^^= , •'
1
_____ U(8,t)
1

CONTINUATION REGION
I TERMINATION REGION
r*(8)
SOLUTION AT TANGENCY
V(8,r)FOR FIXED 8
CONTINUATION
REGION
 TERMINATION REGION
1
V = U ________
V(8,r> ^_____ :5 ==>^'
1
^^^Ui6,r)
1
1
1
SOLUTION AT DISCONTINUITY
OF DERIVATIVE
V(8,t)F0R FIXED 8
CONTINUATION REGION
TERMINATION REGION
V(8,r) ^ — 
V=U
U(8,r)0
SOLUTION AT DISCONTINUITY OF U(8,r>
Figure 3. Possible types of solution
174 J. M. COZZOLINO
of its solution can then be evaluated using the above boundary conditions. This will be illustrated later
for two special cases. Figure 3 illustrates the tangency condition and the two types of discontinuity
conditions.
The method of optimization discussed here is a continuous stage dynamic program. Other mathe
matically similar problems in the literature, such as that of age replacement policy for stochastically
failing devices which wear out, have been formulated as a dynamic program using a "renewal equation"
form. For example, [11], [15], [16], and [17] use that form. Either form could be used here although this
differential form is thought to have some computational advantage for this problem.
The Distribution of Remaining Test Time
When analyzing the feasibility of burnin testing it is useful to know certain statistics of the unit
test. Perhaps the most important of these are the moments of the test duration which will be con
sidered here. Let t be the remaining test time when the test process is in state (0, t), D(t/6, t) its
density, and D(S/6, t) be the Laplace transform of the density. An equation for D{Sld, t), similar
to the value function equation, can be derived using similar argument. By considering a small time
interval, h, and the events possible within it, the following equation is derived:
D(Sld, t) = e™ D(Sld H , t r ) ixti(rld)h + e sh D(SI6, t) [1 /lh>(t/0)A] for t<t*(9).
By allowing h to approach zero, the limiting form obtained is
dD(Sie,T) _^ {T/e) [D(Sldr) _ DiS i eRTH)]=SDiS/dT) for T<T * ( 0).
The associated boundary condition is
D(S/6,t*) = 1.
As a consequence, the Arth moment of remaining test time, defined as
i
E k (d,T) = \ t k D(t/d, r)dt for it = 0, 1, 2, . . .,
obeys the following equation:
dE k (0, r)
/ao(t/0) [E k (d, t) E k (0 H ,T R )] =kE k .t(8, t)
dT
for A=l, 2, . . ., and t 1 t*(0). The boundary condition is
E k {6, t*)=0 for A=l, 2, 3
The solution for these moments requires the knowledge of t*(0); then, since the (k+\) st equation
involves the Ath moment, the equations for the moments can be solved successively starting with
A=l.
5. EXTENSIONS OF THE FORMULATION
The Optimal Policy Under Discounting
In the analysis of processes which take place over an appreciably long time, thus tying up capital
equipment and possibly delaying cash flows, it may be useful to account explicitly for "value of time"
by optimizing the present discounted value of the process. Under this method, the values of all future
OPTIMAL BURNIN TESTING ' 175
value transactions are discounted from face value according to the length of time into the future at
which they take place. Hence, if the amount, £/, is to be received at time, T, into the future, its present
discounted value is Ue~®'\ where /3 is a positive discount factor.
The optimal policy and value function under discounting will be different from that already dis
cussed; however, only a slight modification of the optimization equation is required to treat this case.
The equation now appears as
0= Max {U(d, t)  V(6, T), dV( ° ,T) pV(6, r)  u (tI6) [V(0, t)  V(0 H , r^C,] + C 2 }
OT
The one additional term is (3V(6, t). All that has been said about the nondiscounted optimal solution
applies to this case also.
The effect of present value discounting will usually be to shorten the test. This can be expected
since the terminal value is always received at the end of the test, whereas the test costs are spread
out over the duration of the test. Therefore the expected present discounted value of a given policy is
necessarily less than the expected undiscounted value of the same policy.
BurnIn Testing Under Uncertainty of the Failure Process
We have so far considered a known failure process. However, in many applications there exists
a high degree of uncertainty about the underlying failure process. One method of relaxation of this
assumption of complete knowledge of the density is to consider a known form of density function but
with some unknown parameters. The state of knowledge about the values of such parameters is ex
pressed, at any given time, by a joint density function over all allowable sets of their values. This prior
density is modified by the observation of results of the testing process according to Bayes' Rule. The
concept of a conjugate prior density, as described by Raiffa and Schlaiffer [25], is used here to enable
a mathematically tractable formulation of this problem.
We suppose that the parameterization of the density includes parameter(s) cf> which are unknown,
but do not change with failure. They are now included explicitly in the density, f(t/6, t, </)). The "prior"
density, that which describes the state of knowledge at a given time, is denoted g((/)/i//) . It is normalized
over the range of possible <f) values for any values of its parameters, t/.
The density of fife, unconditional upon the unknown 0, is
Jail <*
h(tie, r, <//) = f(t/e, r, 4>)g{<t>l<\>)d<i>.
Jail <b
The prior density is chosen to be conjugate to f(t/<f), t, </>) so that, after observation of a failure at
time, t , following startup of a unit of effective age, t, the posterior density, according to Bayes' Rule.
/(tie, r, c/QgW)
h(tld, t, i//)
has the same functional form in as does g{<j), <//), but with parameters changed to «//f = t///. (6, t, (//, t).
Hence the value of the parameters of the prior density is sufficient to describe the state of knowledge
at any time.
The failure rate function which appears in the equation is
V(tld, T, (//)=
/;
h(x/6, t, ty)dx
176 J M. COZZOLINO
The state of the process under uncertainty is given by (6, t, i//, x) , where ^ andT are the values current
at the time of last failure or startup and x represents time since that event. Let this description be
abbreviated by the definition
d' = {6, t, ¥}.
The state after failure in state (8' , x) is
d' R {6\ x) = {6, t + x), Tr{6, t + x), V F {d, T, ¥, x)}
x R (d', x) = 0.
It is now apparent that the adaptive burnin test problem, for one unit test, has the same form of opti
mization equation as previously described:
0=Maxj £/(0\ x)V(0', x), dV{ ®' ,x) v(xld')[V(0', x)V(d' R , 0) Ci]+&] .
It will, however, be more difficult to solve because of the larger state space.
6. SOME SPECIAL CASES
There are some important special cases of burnin test processes for which there are no 6 type
variables present; i.e., the failure density of a repaired unit is of the same form as that of some aged
unit which never failed. In such cases, the effective age of the unit (t) completely describes the satte
of the process and the differential equation holds on a line. We now consider three such cases which
encompass a range of possible repair rule cases.
Complete Replacement or Unrepairable Equipment (t« = 0)
In the case of unrepairable equipment, a failure necessitates a complete replacement of tlv* .m J
device by a new, untested one. Therefore the state of the process after failure is independent of the
age attained by the old unit just prior to failure. This case is represented by the repair rule Tr(t) = 0.
The resulting linear first order differential equation,
^^ (jlo(tI6)[V(t) V(0)]^ fil{Tl6)C 1 + C 2
is easily solved to give a family of all solutions parametric on one constant. For convenience, let
w(t) =j „ , ■ +
C 2 f [lFo(xld)]dx
Jo
[lFo(Tld)] [lF (Tld)]
Since the constant of solution is determined by the boundary conditions we write it as t*. The solution
is then
( V{t*)  W(t*) + W{r) for ^ r ^ t*
V(t)=\
[U(t) fort*<r.
The value of t* is then determined as indicated by the previous theorem.
The independence of interfailure times allows easy solution for the Laplace transform of the density
of remaining test time. It is
D{S,t) = Q(s',t*) ''*'*'* f ° r T ^ T *
OPTIMAL BURNIN TESTING J 77
where
1J e Sx Mx/e)dx
Q^,r)=— [ [_ Fo(Tle)]
Although the solution for this case could easily be obtained by more direct arguments, it is pre
sented here both to illustrate the method and to enable a comparison between unrepairable and re
pairable devices. We can already see that burnin testing is somewhat inefficient for unrepairable
devices since a failure results in the loss of all test effort expended up to the time of failure.
Superficial Failure (tr — t)
In the previous case, the test time "invested" in the current unit was completely lost with a failure.
It represents one extreme case the opposite of which is called here superficial failure and results
in no loss of invested time. This case, represented by the repair rule t«(t) =t, is important because of
the simplification of the decision structure.
The differential equation for this case is
The optimal value function in terms of t*, is
\U{t*)C, \ fXo(x/d)dxC 2 (TT) forr^r*
V(t)=\ h
[U(t) for t s*t*.
The Laplace Transform for the density of remaining test time is easily found to be
Z)(S/T)=e s(T *"forTS=T*.
This is the transform of a unit impulse function located at t = t; consequently the length of remaining
test time is not random. The sequential decision which remains unchanged in the face of partial test
experience.
At first glance it seems counterintuitive that a failure leaves the device no better or worse than if
no failure occurred; however, there are reasonable assumptions which lead to this repair rule.
Kooharian and Young [28] present a class of such models. Similarly, Cozzolino [7] discusses a special
case of that class, based upon the concept of production defects initially present in the device. Bazovsky
[4] gives a similar debugging model. Because of its special deterministic nature, the decision policy
is easy to implement. Also, as we shall see, this model results in a relatively profitable burnin test.
Partial Renewal (tr = (xt).
The repair rule tr(t) =oct for 0<a<l represents an intermediate degree of replacement.
Consider, for example, a device composed of N indistinguishable and essential components. Failure
results from failure of any one of the components. Repair is then accomplished by replacement of that
failed component by a new, untested component. The test time previously accumulated on all other
components is retained. This situation may be represented by a partial loss of effective age, where
a in the repair rule is a fraction of the time retained; I J in this case. A more detailed discus
sion of this model is given in [7].
178
J. M. COZZOLINO
The differential equation for this class of repair rules is
dV(i
£ L VL (Tl8)[V(T)V(aT)]=Vo(Tld)C 1 +C 2 .
The combination of differential and difference operators present makes analytic solution difficult for
most failure rate functions of interest. However, approximate solutions based upon numerical integra
tion are easily obtained. Since all solutions are within an additive constant of each other, the boundary
conditions are easily satisfied.
A Specific Example
Some numerical results, based upon the above family of repair rules, were obtained in order to
illustrate how the optimal test policy depends upon the repair rule parameter, a. The specific model
used is based upon the density function,
 b t a \ 6+1
Mtla, b)= 3 *2=0, 6>0, a>0.
The failure rate associated with this density is monotonically decreasing,
b
/xo(tla, b) =
a + t
The value of a tested unit is taken to be a function, R(t), of its first life after test, t, as shown in Figure 4.
I 15 Ro'
Ro'
85 Ro
RO) ' VALUE RECEIVED FOR TESTED UNIT
NORMAL
INCENTIVE RANGE
PENALTY
RANGE
RANGE
Figure 4. Value vs performance
The reliability value function is consequently the expected value of this function as of test termination
and is a continuous function of age,
U(t)
\:
R(t)f(tl6,r)dt
= R«
0.85 + 0.15
a + T
a + T + 7 7 ,
+ 0.15
a + r
a + T + T 2
The parameter values chosen for the example are a = 72.0, 6=1.1, Ci = 1.0, C 2 = 0.01, R o = 90.0,
ri = 500, and T 2 = 1000. The mean life is found to be 720 units and the value of an untested unit is
U(0) = 78.57. The optimal test policy and value function were found for several values of repair rule
parameter, a. The optimal policy, r*, is tabulated for these a values in Table 1, along with other quan
tities of interest. Figure 5 gives the optimal value function for one a value; a = 7/8. Figure 6 shows how
the expected test profit, V(0) —U{0) increases with a.
OPTIMAL BURNIN TESTING
179
89

87

^
Sale Price
85
83
>Test Cost
81
79
/
■
>Tesl Profit
77
I 1 1 1 1 1 1 1
1 1 1 1 1 1
40 80 120 160 200 240 280 320 360 400 440 480
FIGURE 5. The optimal value function V{t) for a = 7/8
V(0)U(0)
18
16
14
12
1.0
08
06
04
02
n
 t — i — i— ~~\ i i i i i
c
) 01 02 03 04 05 06 0.7 0.8 09 1
D
6 8 10 =
la
i iii mi i
Figure 6. Test profit as a function of a;
Case 3
The parameter, a, indexes a family of partial renewal models; the fraction of the device which is
renewed with repair of a failure is (1 — a)= 1/iV. As a ranges from to 1 the model ranges from com
plete replacement (a=0) to the superficial failure case a— 1, and is therefore capable of representing
a wide range of repair behavior. The computations show that as the degree of partial renewal decreases
{a increases), the coefficient of variability of the total test duration decreases, as would be expected
from previous results. As the test duration becomes more deterministic the optimal test policy, t*,
lengthens. The expected profit due to testing increases and, in addition, the ratio of expected test
profit to expected test cost increases, showing an inherently greater profitability as the renewed portion
of the device decreases. This result conforms to our intuitive expectation based upon the previous
results; as a increases a larger portion of the invested test time is retained at failure. Consequently,
failure becomes less of a risk and the test process becomes closer to a deterministic one.
180
J. M. COZZOLINO
Table 1. Optimal Test Policy for a Class of Repair Models
a
N=aar>
Optimal
policy
T*
Expected value
Expected
profit from
test
[V(0)U(0)]
Ratio (%) of
test profit
to test cost
r V(0)U(0) i
IV(T*)V{0)\
Expected
duration of
unit test
Coefficient
of variability
of duration of
unit test
Before test
V{0)
After test
V(T*)
1
14.5
78.595
78.979
0.0249
6.5
15.7
0.28
0.5
2
80.0
78.70
80.68
0.13
6.6
97.5
0.26
0.6
2.5
140.0
78.83
82.04
0.26
8.1
174.2
0.25
0.7
3.3
220.0
79.09
83.63
0.52
11.5
268.4
0.20
0.75
4
261.0
79.27
84.36
0.70
13.8
310.9
0.17
0.8
5
304.2
79.47
85.07
0.90
16.1
351.7
0.14
0.875
8
370.0
79.82
86.05
1.25
20.1
406.4
0.09
0.9
10
390.0
79.94
86.32
1.37
21.5
420.6
0.07
1.0
oc
470.3
80.43
87.36
1.86
26.9
470.3
0.0
7. CONCLUSIONS
The analysis of the burnin testing of repairable devices requires explicit information about the
repair process in order to evaluate repaired devices. A class of repair models, having a form tractable
with respect to the sequential decision problem for optimization of the burnin test, is described by a
repair rule of the form, 0«(0, t), tr(0, t), specifying the parameters of a density function within the
family of the renormalized form of the new unit density. This family of densities has the same func
tional form of failure rate.
The burnin testing optimization, when treated as a continuous stage dynamic program, results
in a differential equation similar to equations for the state probability of a stochastic process. In two
important special cases, the differential equation is of the first order in total derivatives and conse
quently is easily solved and fit to the boundary conditions. In more complex cases it has a differential
difference nature.
The wide differences in the expected test profit under different repair models of the partial renewal
family leads to the conclusion that the explicit modeling of the failurerepair process is highly important
for the burnin testing of repairable equipment. While this result is based upon a specific class of
models, it appears likely that the greater profitability of burnin testing of repairable over unrepairable
devices may prevail over a much wider class of models.
REFERENCES
[1] Barlow, R. E., A. W. Marshall and F. Proschan, "Properties of Probability Distributions with Monotone Hazard Rate," Annals
of Math. Stat., 34, 375389 (1963).
[2] Barlow, R. E. and F. Proschan, Mathematical Theory of Reliability (Wiley, New York, 1965).
OPTIMAL BURNIN TESTING  181
[3 Barlow, R. E. and E. M. Scheuer, "Reliability Growth During a Development Testing Program," Rand Memorandum
RM4317NASA (Nov. 1964).
[4] Bazovsky. L, Reliability Theory and Practice (Prentice Hall Inc., New Jersey, 1961).
[5] Chemical Rubber Publishing Company, Standard Mathematical Tables, 12th edition. Commissioner's 1941 Standard Ordi
nary Mortality Table.
[6] Cox, D. R.. Renewal Theory (Methuen and Company, Ltd.. London, 1962).
[7] Cozzolino. J. M., "The Optimal BurnIn Testing of Repairable Equipment." Massachusetts Institute of Technology, Opera
tions Research Center Technical Report No. 23 (Oct. 1966).
[8] Davis, D. J., "The Analysis of Some Failure Rate Data," J. Am. Stat. Assoc, 47, 113150 (1952).
[9] Drenick, R. E., "The Failure Law of Complex Equipment," J. Soc. Indust. Appl. Math., Vol. 8, No. 4 (Dec. 1960).
[10] Dreyfus, S.. "Dynamic Programming and the Calculus of Variations," Rand Dept. R441PR.
[11] Fox, B. L., "An Adaptive Age Replacement Policy." Operations Research Center, University of California, Berkeley 6517
(PR) (1965).
[12] General Electric Co., Transistor Manual, Chap. 18, "Semiconductor Reliability," 6th Ed. (1962).
[13] Haddon, M. O, H. Cheesebrough, and R. E. Kirby, "Guaranteed Reliability at a Minimum Cost," IEEE Trans, on Reliability
and Quality Control (Sept. 1958).
[14] Hausman, W. H., and M. Kamins, "Early Failures in Automobile Parts: A Background Study in Reliability," Rand Mem.
RM4002PR (1964).
[15] Jorgenson, D. W. and J. J. McCall, "Optimum Scheduling of Replacement and Inspection," Operations Research II, 732746
(Oct. 1963).
[16] Jorgenson, D. W., and R. Radner, "Optimum Replacement and Inspection of Stochastically Failing Equipment," in Studies
in Applied Probability and Management Science (Stanford University Press, Sanford, 1962), edited by K. J. Arrow, S. Karlin,
and H. Scarf.
[17] Kamins, M. and J. J. McCall, "Rules for Planned Replacement of Aircraft and Missile Parts," Rand Mem. RM2810PR
(Oct. 1961).
[18] Knight, W. R., "Exponential and Subsexponential Distributions in Statistical Life Testing," Ph. D. Thesis, University of
Toronto (1958).
[19] Lawrence, M. J., "An Investigation of the BurnIn Problem," Technometrics, 8, 6171 (Feb. 1966).
[20] Lomax, K. S., "Business Failures: Another Example of the Analysis of Failure Data," J. Am. Stat. Assoc, 49, 847852 (1954).
[21] Markowitz, D. H. and R. G. Pouliot, "Infant Mortality Screening; A Discussion of BurnIn and Other Screening Techniques,"
U.S. Naval Weapons Quality Assurance Office, Washington, D.C. (June 1967).
[22] Morse, P. M., Queues, Inventory, and Maintenance (Wiley, New York, 1958).
[23] Norris. R. H., "RunIn or 'BurnIn' of Electronic Parts," Proceedings of the Ninth National Symposium on Reliability and
Quality Control (1963), pp. 335357.
[24] Proschan, F., "Theoretical Explanation of Observed Decreasing Failure Rate," Technometrics, Vol. 5, No. 3 (Aug. 1963).
[25] Raiffa, H. and R. Schlaifer, Applied Statistical Decision Theory (Harvard University Press, Cambridge, 1961).
[26] Riordan, J., Stochastic Service Systems (Wiley, New York, 1962).
[27] Watson, G. S.. and W. T. Wells, "On the Possibility of Improving the Mean Useful Life of Items by Eliminating those with
Short Lives," Technometrics, 3, 281298 (May 1961).
[28] Young. H. and A. Kooharian, "A Decision Model for the Operational Acceptance Test (OAT) of Minuteman Ground Elec
tronic Equipment," Sylvania Electronic Systems, Res. Rept. No. 477 (Sept. 1965).
[29] Zelen. M. (Ed.). Statistical Theory of Reliability, University of Wisconsin Press, Madison, 1963.
GOAL PROGRAMMING IN ECONOMETRICS
W. Allen Spivey
University of Michigan
and
Hirokuni Tamura
University of Washington
1. INTRODUCTION
Optimization models have been suggested as an alternative to simultaneous equation models for
use in dealing with the problems of policy formulation and decision making at the level of the national
economy (Theil [11]; Holt [5]; van Eijk and Sandee [14]; and Hickman [4]). In simultaneous equation
models a policy maker prescribes a set of fixed values for goal or target variables and these are inter
related with a set of instrument variables by means of a system of simultaneous linear equations. The
policy maker then solves this system (in which the number of equations or restraints is usually required
to be equal to the number of unknowns) for the values of the instrument variables. It has been empha
sized (Tinbergen [12], [13]) that the policy maker's welfare function is not explicit in this model but is
implicit in the selection of the goal and instrument variables. In the optimization approach one develops
a welfare (preference) function which is optimized with respect to a set of econometric restraints. As
is usual with such models, there can be many feasible solutions, and the optimization process then
chooses a best or optimal solution from the set of feasible solutions.
In the Theil optimization model (Theil [11]) a preference function is developed in which both goal
and instrument variables appear. A desirable configuration of goal and instrument variables, not neces
sarily satisfying the constraints, is specified by the policy maker, and a quadratic function in deviations
between the desired and feasible values of the variables is optimized (minimized). Although the quad
ratic function allows the policy maker to take account of a decreasing marginal rate of substitution
between any two variables (at least over an appropriate subset of the domain of the function), it has the
disadvantage of assigning the same penalty to overattainment of a goal as to underattainment (Hick
man [4], p. 6). A policy maker, for example, might not want to accept underattainment of a goal, but
might be indifferent as to whether the goal was met or exceeded.!
In this paper, we introduce an optimization model which differs from that of Theil in several ways.
It permits overattainment and underattainment of any goal, yf, to be weighted either equally or dif
ferently in the decision maker's preference function. The model can be regarded as a special type of
linear optimization model (as we make clear below) so it is easy to handle computationally and a wide
range of parametric or sensitivity analyses can be implemented on a computer. These analyses can
tStill another approach has been developed (Spivey and Tamura [9]), that of using a generalized simultaneous econometric
model (GSEM), where the word generalized is used to emphasize that the number of equations or constraints is not required to
be equal to the number of variables. The GSEM can be regarded as a link between conventional simultaneous equation models
and the Theil model since, like the former it consists of simultaneous linear equations and like the latter a GSEM can be shown
to have a solution which is identical with an optimal solution to a Theil type optimization model.
183
184 W. A. SPIVEY AND H. TAMURA
display interactions between alternate optimal policies and goal attainment and can assist the policy
maker in studying explicitly many of the trade offs between underattainment, attainment, and over
attainment of selected goals. Such sensitivity analyses are not available for either the Theil model
or Tinbergen type simultaneous equation models. The model can also easily accommodate explicit
restrictions on instrument variables in the form of nonnegativity and upper and lower bounds.
The model is that of goal programming, originally developed by Charnes and Cooper in another
context (see Charnes and Cooper [1]; also Charnes, Cooper, and Ijiri [2]; Ijiri [6]). We develop an appli
cation of this model using Klein's Model I of the United States economy as a basis (Klein [7]) and display
in numerical form and interpret a sequence of decision making implications of several classes of the
sensitivity analyses that can be performed with this model.
2. GOAL PROGRAMMING AND GOAL ATTAINMENT
A reduced form of a linear econometric model can be written
(1) y=Rx + s,
where y is an n by 1 vector of goal or target variables, x is an m by 1 vector of instruments,/? is ann by m
(real) matrix of impact multipliers, and 5 is an n by 1 vector of constant or additive terms. If m — n, this
is the reduced form of a conventional system of simultaneous linear equations; if m ^ n, then (1) is an
example of a generalized simultaneous econometric model, in which m = n is a special case (Spivey and
Tamura [9]).
Let y* be the goal vector prescribed by the policy maker and consider the model
(2) y* = /?* + s,
where R and 5 are assumed to be known. If R is singular, then (2) either has no solution or infinitely
many solutions. If R is nonsingular, then (2) always has a unique solution; however, it may not be
acceptable to the policy maker because one or more components violate restrictions which are present
in the mind of the policy maker, but are not explicitly a part of (2) — for example, nonnegativity.t
We reformulate (2) as
(3) y*s = RxIz + + h
x,z + ,z~mO,
where / denotes an n by n identity matrix, R is the n by m matrix introduced in (2), and z + and z~ are
n by 1 (unknown) deviation or discrepancy vectors. The ith components zf of z + and zf of z" measure
the deviation upwards and downwards, respectively, of feasible y, values from the corresponding goal
values yf; in other words, the ith constraint of (3) can be written
(4) yts^RfXzt + zi (»=1, ...,»),
where /?, denotes the ith row vector of the matrix/? and s, the ith component of the vectors. Note that
a nonnegativity restraint on x is explicit in this formulation.
If the policy maker can find a feasible instrument vector x for which bothz+ andz,~ in (4) are zero,
then the goal yf can be attained exactly. On the other hand, if there is no feasible vector x for which
zt—zf = in (4) he cannot attain the prespecified goal yf. He can, however, find a feasible x that will
tlf implicit restrictions on x are violated in the case R is nonsingular, Tinbergen [12], [13] would suggest model reformulation
by selecting another goal vector y*' in an ad hoc manner, substituting y*' for y* in (2). and again solving for x.
GOAL PROGRAMMING 185
allow him to come as close as possible to y* in a sense we will soon make clear. Now consider the
following goal programming model
(5) minimize e T z + + e T z~,
subject to
(6) y*s = Rxh + + Iz
x,z + ,z^0,
where e T is an n by 1 row vector t each of whose components is 1.
The z+ and zf (i= 1, . . . , n) enter into the objective function (5) with a positive coefficient so the
minimization process reduces the values of these variables in an optimal solution to the lowest feasible
level and (5) attains a minimum value of zero if and only if both z + and z~ are the zero vectors (if and
only if zt=Zi~ = Q for each i). The equations (6) clearly have feasible solutions because the rank of
the matrix / assures that (6) can be solved for the z's alone.
The advantage of reformulating (2) as a goal programming problem should now be clear: if (5)
has a minimum value of zero, then (2) has a solution vector x which enables the policy maker to attain
his prespecified goals yf (1 = 1, . . ., n) exactly. If it does not, then (2) is not solvable for a feasible
instrument vector x. In the latter case, however, an optimal solution to (5) and (6) will disclose (a)
how close the policy maker can get to the desired goal vector y*, and (b) whether this approximation
is on the high side or low side for each individual goal, y* .
The statement (b) can be made, and nonsense values of zf and zj avoided, because the vectors
in the system of restraints (6) associated with zf and zj are negatives of each other. It is impossible
for both zf and z,~ to assume positive values in an optimal solution because the basis property of
linear independence is preserved in the simplex solution procedure. Denote an optimal solution to
(5) and (6) as
'Xo
(7)
where xj= (*i,o v • ., x m ,o), Zq T = (zt, o> • • •'> z n, o) an( ^ z o T= ( z i~, o» • • •> z n , o) There are exactly three
possibilities in an optimal solution to (5) and (6):
(i) zt,o = Zi,o — 0,
(ii) <•> 0,^0=0,
( iH ) z<to=0, zf, >0.
t The superscript T denotes a transpose. All vectors will be regarded as column vectors unless this superscript is used. The
model (5) and (6) can also be written
Minimize [0, e T , e r ] z*
Subject to V
y*s= [/?,/, /] lz<
x
z* \
186 W. A. SPIVEY AND H. TAMURA
Case (i) means that the ith goal y* is attained exactly whereas cases (ii) and (hi) indicate that a deviation,
upwards and downwards respectively, must be accepted by the policy maker desiring a nonnegative
solution to (2).
The loose phrase, "how close the policy maker can get to the goal vector y*," is given a precise
meaning in goal programming because in (5) and (6) one is minimizing
W X \yf(Rix+ Si )\,
i
where the vertical bars denote absolute value. An optimal solution x yields values of the discrepancy
variables Z+ and zj which have the following effects:
(i') Zi,o = z7,o =Q if y* = RiXo + sr,
(ii') zt >0 =RiXo+Si — yf, z f " =0 if y?<R(Xo+Si;
(iii') Zi,o=y*— (Rtxo + Si), zl o =0 i( y*> RiXo + si.
This can be stated in still another way. Given an xo in an optimal solution (7), the corresponding target
configuration y , which occurs when xo is used in (2), is
(9) y = Rx + s.
If case (i) occurs, yo = y* as desired, but in either case (ii) or (iii) y ¥= y*. In the latter cases yo is the
closest vector to y* where closeness is measured by the absolute value metric or distance measure
(this metric is sometimes called the L t metric). That is, yo is an optimal solution to the problem
n
(10) minimize^ \y*~ Yi\,
subject to
(11) y = Rx + s,x^0.
The measure of closeness can be generalized by using an L p metric and minimizing
(12) (2ly<** p V
subject to the restraints (11), where p is any integer greater than or equal to one. If p = 2, this is the
same as minimizing the quadratic function X(y* — y;) 2 subject to the restraints (11), which is the form
of the objective function in the Theil optimization model. The model in this case selects a feasible
vector y which is closest to y* in the sense of Euclidean distance or ordinary least squares.
The goal programming model (5) and (6) can be augmented so as to contain political and institu
tional as well as nonnegativity restrictions on the instrument variables x,(i = l, . . ., m). Such
restrictions, called boundary conditions by Tinbergen ([12], [13]), must be dealt with implicitly in
the use of his models, but can be incorporated as upper and lower bounds on xi. These offer no difficulty
in a linear programming model. For example, we could
(5) minimize e T z + 4 e T z~
subject to
(6') y*  5 = R X  Iz + + Iz
GOAL PROGRAMMING 187
Xi^ Li
US,
: k ^ U k
keT,
X, z + , z~
i=0,
where S and T are any subsets of the set of integers {1, . . ., m} and Li and Uk are specified lower
and upper bounds, respectively. The model (5) and (6') can be regarded as a decision model which
explicitly takes account of all relevant restraints, both econometric and political or institutional.
3. WEIGHTING DISCREPANCY VARIABLES
In addition to determining a feasible vector y which is closest to the desired vector y* in terms
of the L\ metric, a policy maker might regard attainment of certain goals as more important than that
of others and might wish to use different weighting factors on deviation variables to reflect this. One
can accordingly replace e T in the objective function (3) by other nonnegative weighting vectors u T
and v T . This gives the following objective function instead of (5),
(13) u T z + + v T z~
which is minimized subject to (6). Also,
u T = v T = e T
is a special case in which each deviation variable is given a weight of one. Thus if a policy maker
desires to attain a given goal yf more than others he can attach larger weights m and v ,• to the corre
sponding deviations zf and z,~, respectively. If he would accept an overattainment of yf, but cannot
tolerate any underattainment, he could set Uj = 0, but would make Vj large enough to prevent zj from
assuming a positive value in an optimal solution. t In this case, the model will allow the policy maker
to impose the following condition
(14) yf =i R jX + s h
i.e., we can accommodate an inequality restraint if desired — a feature which is not incorporated in
the Theil optimization model.
Moreover, one can parameterize the objective function in a goal programming model to see how
an initial optimal deviation package — say, one in which some goals are attained exactly and some are
not attained initially — varies as different nonnegative weights are attached to the deviation variables.
This can be contrasted with the effect of a weighting factor in a Theil model (L 2 metric model), which
permits a delicate control of the size of the corresponding deviation symmetrically above or below
the goal vector y* (for an empirical study of the difference in effects of the use of Li metric and L 2
metric, see Tamura [10]).
4. SOME EXAMPLES OF DECISIONS IN GOAL PROGRAMMING MODELS
The classical Klein Model I (Klein [7]) will be used to present some numerical examples of goal
programming decisions. Model I in reduced form is shown below; actual values have been given to the
lagged variables and the trend variable so that model can be used in examining 1933 policy decisions. tt
tFor any nonnegative vector of coefficients [u T v T ], (13) has a lower bound of zero in a minimization problem. Only models
with nonnegative coefficients in the objective function will be considered in this paper.
ttThis model was also used for illustrative purposes by Theil [11] and Spivey and Tamura [9].
188
Consider
(15)
W. A. SPJVEY AND H. TAMURA
0.6659
0.1883
0.6705'
0.2238
1.2808
1.1188
0.1620
0.2033
0.8102
0.0517
0.2959
0.2584
0.6142
1.4842
1.9289
0.0517
0.2959
0.2584
+
/37.58^
' 9.0764
23.0109
5.4938
32.0873
i 201.6060/
where C= consumption expenditures, P= corporate profits, W\ — private wages, /= private invest
ment, Y— national income, and K= end of year capital stock are endogeneous variables and W<l = gov
ernment wages, 7"= tax revenue, and G = government expenditures are exogeneous variables (all
variables are measured in terms of billions of 1934 dollars). The constant terms are calculated under
the assumption that the random disturbances are forecast accurately.
EXAMPLE 1. Parameterization of Weights: Suppose that the goal vector selected by a policy
maker is y T = (C, P, Wi, I) and the instrument vector is x T = {W2, T, G,). Then in the manner of Tin
bergen and Theil we can formulate a linear econometric model in these vectors y and x by extracting
the first four equations from the reduced form (15).
Let the goal values be set as a uniform 10 percent increase over the actual values of 1932 for each
goal variable, i.e., y* T — (50.16, 7.7, 31.9, —5.58). Then we use the following goal programming model
to examine 1933 policy decisions,
(16)
subject to
(17)
minimize zt + z c + zp~ + z p + z^ x + z Wl + zf + Zj
+
0.6659X
0.2238 \ w
0.1620 I
0.0517/
i\ , / 1
0/
0.6705\
1.1188 G +
0.8102
0.2584
U*i +
W 2 , l , tr, Z c , Z c , Z p , Z p , Z Wl , Z Wi , Z/ , 2/ = U,
where (16) and (17) correspond to (5) and (6), and where each goal is given equal weight. An optimal
solution to this problem leads to the policy
(18)
6.870 1
I 13.650,
,15.770/
and the corresponding target configuration yo is
(19)
yo
GOAL PROGRAMMING
189
Comparison of yo with the prescribed vector y* indicates that the three goals C*, P*, and W* are
attained exactly and the private investment goal /* is underattained by $.23 billion.
Suppose that the policy maker, upon being informed of (19), wants to assign a larger weight to the
attainment of the investment goal/* than is the case in the objective (16). Then the deviation variable
zj can be given the weight A (z/ being given the same weight 1), and we
(20)
minimize z£ + z c + zp~ + z p + z+ 1 + z Wi + z/" + kz,
subject to (17). We study this optimization for various values of k (continuous parameterization of
k). Table 1 summarizes relevant information. For 1 Si k < 4.328 an optimal solution remains unchanged
and for each such choice of weighting factor k, Xo in (18) is optimal and y in (19) is realized. This is indi
cated by the first data row of Table 1. For weight k ^ 4.328, we must shift to a new optimal solution x ,
getting a corresponding y . The second data row in Table 1 displays values of these vectors. The invest
ment goal is attained for weights k in the range indicated and we overattain the goal P* by one billion.
Table 1. Optimal Policies and Consequent Goal Attainments
Range of Values
of A
Components of y
Components of xo
Co
Po
r,,
h
w 2 ,
T
Co
lgK 4.328
50.16
LI
31.9
5.71
6.870
13.650
15.770
4.328 S k < °°
50.16
8.7
31.9
5.48
6.845
12.641
15.512
Note: Underlining indicates exact attainment of goal.
Comparison of the two optimal policies indicates that attainment of /* is brought about by in
creased deficit spending; Go  To = 15.770— 13.650 = 2.120 billion for the original optimal x , whereas
G — To = 15. 512— 12.641 = 2.871 for the alternate optimal xo Note that the increased deficit is not
attributable to increased values of Go, but rather to decreased values of T . That is, in this model
optimal adjustments call for reducing taxes rather than increasing government expenditures.
A policy maker might wish to inform himself further about the interactions between T, G, I, and
optimal policies. T, G, and / are related by means of the fourth equation in (15),
/ =  0.0517r 2  0.2959r+ 0.2584G  5.4938.
One could study how sensitive an optimal policy is to changes in either of the multipliers — 0.2959 and
0.2584. A computer program can be used to parameterize the matrix in (17), one row or one column at
a time, so parameterizing any one element in the matrix is a special case. Specifically, one could
determine a range of values containing — 0.2959 within which an xo remains optimal, assuming that
other columns of the matrix are held constant. Similar remarks can be made concerning the para
meterization of the government expenditures multiplier 0.2584. Indeed, the following kinds of para
metric analysis can be performed on the model (16) and (17): (a) nonnegative parameterizing of the
coefficients in (16); (b) parameterizing the column vector of constants in (17); (c) parameterizations of
(a) and (b) simultaneously; (d) parameterizing the matrix in (17) one row or one column at a time. Thus
the model provides the policy maker with considerable flexibility in both policy formulation and decision
making.
190 W. A. SPIVEY AND H. TAMURA
EXAMPLE 2: Prespecification of Both Goals and Policies: Theil [11] introduced an optimization
model in which the policy maker specifies a goal vector, y*, and a desired instrument vector, x*. We
reformulate this model by means of goal programming and introduce additional flexibilities by using
alternate objective functions.
The target or goal variables chosen by the policy maker are C, I, and D, where D is a distribution
variable defined as W x — 2P (see Theil [11], p. 79). Prespecified goal values are, in billions of 1933
dollars, C* = 49.69, /* = — 3.10, and D* = 11.25. The desired instrument values for W 2 , T, and G are
5.04, 7.40, and 10.44, respectively (for specification of goal and instrument variables, see Theil [11]).
Thus
(21) y* T = (C*,I*,D*) = (49.69 3.10 11.25)
and
(22) x* T = {Wt,T*,G*) = (5.04 7.40 10.44).
The goal programming model is to
(23) minimize e T z + + e T z~
subject to
(24) /y*\ /s\ (R
X, Z + , 2^0,
, (A .
where I in (24) comes from (21) and (22), and where
/ 37.5812 \ / 0.6659 0.1833 0.6705'
s = 5.4939 I, R= 0.0517 0.2959 0.2584
\ 4.8580/ \ 0.2855 2.3583  1.4274,
were calculated from data in Klein's model 1 in the manner suggested by Theil (see Theil [11] p. 87
and Tamura [10], p. 46), h and h denote identity matrices of size 3 and 6, respectively, and
Z + = (Zf , Zj , z D , z W2 , Z T , Z (: ) ,
z~ = {zj}, ZJ, zjj, z W2 , Z T , z a ) ,
x T = (W 2 , T,G).
As observed earlier, the objective function we are using assigns the same weights to each zj 1 " and
zf. If the policy maker would accept an overattainment of goals C* and /* and an underattainment of
goal D*, then minimization of the following objective function subject to (24) could be suggested,
(25) f
where/ is a column vector which contains 12 elements all but 3 of which are set equal to 1. The remain
ing three elements are the coefficients of zf , zf, and zj, and are set equal to zero. The zero coefficient
of these variables then essentially excludes them from the minimization process.
GOAL PROGRAMMING
191
Optimal policy solutions under the two different objectives are given in Table 2. For purposes of
comparison we have also indicated a Theil type optimization decision, i.e., one from a L 2 metric model.
TABLE 2. Alternate Optimal Policy Decisions (Billions of 1934 Dollars)
Values of y
Values of x
C
I
D
w 2
T
G
Desired Values
49.69
3.10
11.25
5.04
7.40
10.44
Optimal Solution
with (23)
46.35
5.55
11.25
5.04
8.42
10.44
Optimal Solution
with (25)
46.54
5.25
8.84
5.04
7.40
10.44
Theil Model
47.77
5.41
10.42
6.44
8.32
11.13
Note: Underlining indicates exact attainment of goal.
5. CONCLUDING REMARKS
Earlier we indicated a generalization of optimization models in terms of the L p metric, a generali
zation which includes both goal programming and Theil type optimization models as special cases, i.e.,
subject to
minimize I ^ \yf — J\  p ) 1/p ,
y = Rx + s, x ^ 0.
Although for p = 2 this is a nonlinear optimization problem and, as such, may appear tQ present
insurmountable computational difficulties, recent research (Duffin, Peterson, and Zener [3] and
Peterson and Ecker [8] ) has shown that this problem is one of a class of nonlinear optimization prob
lems known as geometric programming. Results for these models are available whereby numerical
calculations and approximations can be developed. Moreover, certain types of duality theorems and
"complementary slackness" properties have been established, the economic interpretations of which
are yet to be explored. Research by the authors on these matters is now in progress.
REFERENCES
[1] Charnes, A. and W. W. Cooper, Management Models and Industrial Application of Linear Programming (John Wiley &
Sons, Inc., New York, 1961), Vol. 1.
[2] Charnes, A., W. W. Cooper, and Y. Ijiri, "Breakeven Budgeting and Programming to Goals," J. Accounting Research /, 1043
(1963).
[3] Duffin, R. J., E. L. Peterson, and C. Zener, Geometric Programming— Theory and Application (John Wiley & Sons, Inc.,
New York, 1967).
[4] Hickman, B. G. (editor). Quantitative Planning of Economic Policy (The Brookings Institute, Washington, D.C., 1965),
Chap. 1.
[5] Holt, C. C, "Quantitative Decision Analysis and National Policy: How Can We Bridge the Gap?" in Quantitative Planning
of Economic Policy, edited by B. G. Hickman (The Brookings Institute, Washington, D.C.. 1965).
[6] Ijiri, Y., Management Goals and 'Accounting for Control (North Holland Publishing Company, Amsterdam; Rand McNally
Company, Chicago, 1965).
[7] Klein, L. R., Economic Fluctuations in the United States 19211941, Cowles Commission Monograph No. 11 (John Wiley &
Sons, Inc., New York; Chapman and Hall, Ltd., London, 1950).
192 W. A. SPIVEY AND H. TAMURA
[8] Peterson, E. L., and J. G. Ecker, Geometric Programming: Duality in Quadratic Programming and L v — Approximation,
The University of Michigan, Technical Report No. 19682, 1968.
[9] Spivey, W. A., and H. Tamura, "Generalized Simultaneous Equation Models." to be published in International Economic
Review.
[10] Tamura, H.. "Linear Models for Macroeconomic Policy Making" (unpublished Ph.D. thesis) The University of Michigan,
1967.
[11] Theil, H., Optimal Decision Rules for Government and Industry (Rand McNally and Company, Chicago, 1964).
[12] Tinbergen, J., "On the Theory of Economic Policy," of Contributions to Economic Analysis (North Holland Publishing
Company, Amsterdam, 1952), 2nd ed., Vol. I.
[13] Tinbergen, J., Economic Policy: Principles and Design of Contributions to Economic Analysis (North Holland Publishing
Company, Amsterdam), Vol. XI.
[14] van Eijk, C. J. and J. Sandee, "Quantitative Determination of an Optimum Economic Policy." Econometrica 27, 113 (1957).
PRODUCTION AND EMPLOYMENT SCHEDULING IN MULTISTAGE
PRODUCTION SYSTEMS
Christoph Haehling von Lanzenauer*
The University of Western Ontario
London, Canada
ABSTRACT
A model is developed for planning optimal production and employment levels in multi
product, multistage production systems. The market requirements for each product over
the planning period are given. Backorders and/or shortages are permitted. Backorders
and shortages must be considered in order to determine the amount of each product's demand
that should be filled, backlogged, or go unsatisfied if the production capacity is insufficient
to fill all market requirements. Backorders and shortages, on the other hand, are desirable
under certain dynamic market conditions.
I. INTRODUCTION
The production and employment scheduling problem — planning production and employment levels
over time — has received a great deal of attention. Models and decisionrules generally based on single
stage production systems have been developed for many situations under various conditions [7]. When
ever many products are to be manufactured the concept of aggregate production scheduling has
frequently been suggested [1, 5]. The aggregate production level is then represented by some appro
priate measure of production per unit of time [5] in view of the common denominator required. Once
the aggregate production level is specified the problem of determining the quantities and the timing
for the individual products still remains. This problem of disaggregation is of special importance when
ever the available capacity is insufficient to produce all market requirements and therefore some
demand has to go unsatisfied. It is not at all obvious how the disaggregation decisions should be made
in view of insufficient capacity.
Following along the lines of more recent work in production planning [2, 3, 8], we consider a multi
stage production system consisting of various production facilities (stages) linked together by inprocess
inventories. The objective of this paper is to develop a model for multiproduct, multistage production
systems by which optimal production decisions on a disaggregate basis for both sufficient and insufficient
production capacity can be made. As a point of departure we take Zangwill's [8] model. The subsequent
analysis, however, differs from Zangwill's work on a number of points:
(a) The demand requirements need not be met.
(b) Each production stage manufactures various different products.
(c) Workforce decisions are integrated.
(d) A different cost structure is assumed.
II. THE MODEL
The Production System
The elements of a multistage production system are the individual production stages consisting
of production facilities and an. inventory. These production facilities can be thought of as multipurpose
* Assistant Professor of Business, School of Business Administration.
193
194
C. H. von LANZENAUER
machinery. The output of a given production stage will be fed into the inventory bank from which the
units can be used either as inputs to higher numbered production stages or to supply market require
ments for the stage's products. Figure 1 describes such a production system.
"ff~
MARKET
Production and Inventory Aspects
Let Dj, be the market requirements of production stage fs (j= 1, 2, . . . , J) product i (i= 1, 2,
...,/) during the period t(t=l, 2, . . ., T). It is assumed that the requirements in each period
over the planning period are known in advance. The planning period is medium range so no discounting
is required. The market requirements are expressed in units. Note that the index i at each production
stage represents physically different products. Let Hj, be the amount of product i in the inventory of
stagey at the end of period t. The inventory equation can then be expressed by (1).
(1)
tih = H h
+ X ( +X j
D h
k=j+l n=l
for all i, t andy= 1,
,yi
with XjtR {Xj l0 ) being the amount of product i produced in stage j during period t on regular time
(overtime). aj£ (ajg 3= 0) can be considered as a production coefficient and is the number of units of
production stager's product i required to produce one unit of stage &'s (k=j+l, j+2, . . ., J)
product n (n= 1, 2, . . ., /). aj£ = O if product i of stagey is not being used for the production of
stage Fs product n. The inventory equation for the last stage J of the production system is
(2)
Ui — JJi _l_ Vi _L Yi _ f)i
for all i, t.
The formulation of (1) and (2) assumes that the production capacity available is sufficient to satisfy
all market requirements. If market requirements exceed the available production capacity, the excess
demand must either be backlogged into the next period or has to go unsatisfied. The possibility for
backlogs and shortages is only meaningful for market requirements since the requirements of sub
sequent production stages, i.e., the production levels, are variable themselves. The inventory equation,
adjusted by backlogs and shortages, can then be expressed by (3)
(3) n}, = Hj t _ l +X} at +X} lo (Dl t + B} t _ l )+Bj t +Sl t  ^ z, » *"*m > "tto
k=j+] n=i
for all i , t and y = (/' + 1 ,
J I)
with #j ( being the demand backlogged into period £+ 1 and Sj, the amount of the periods total demand
(Dj t + Z?j,_,) which will never be filled. Note that the variables//],, Bj t , and SL are nonnegative vari
ables. Generally, one can assume that, whenever a positive inventory exists of production stage y's
product i at the end of period t , Bj, and Sj, will be zero; on the other hand, whenever #j, and/or Sj, is
greater than zero, there will normally be no inventory available. It should be realized, however, that
the formulation of (3) permits the possibility of having positive inventory as well as backlogged and/or
unfilled demand at the same time. This configuration can be desirable considering over time changing
MULTISTAGE PRODUCTION SYSTEMS 195
economic conditions. For example, in view of a speculative increase in the selling price in periods
ahead, it might be advantageous to let some of the current periods demands go unfilled, though inventory
is available. The inventory equation for the last production stage is
(4) //', = //),_, + X> m + X* M  (/>*,+ 0j,_, ) +B< ( + Sj, , for all i, t.
For the special case of a production line (i.e., production stager's output will only be used in produc
tion stage;' + 1) we have instead of (3),
(3a) ^ = ^ 1 +^^o(^ + ^ 1 )+^ + ^2 «&, C*£„»+*W'
n=l
for aU i, t.
In view of hmited inprocess inventory capacities, we can restrict
(5) EP/ffJ* 1 * for all;,*,
i=\
where /3j is the amount of space required to store one unit of production stager's product i. Lj t repre
sents the available inventory capacity in space terms. Correspondingly, restrictions can be placed on
the maximum amount to be backlogged per period. Such an upper limit is desirable for Bj T in view
of the fact that the cost associated with backlogs are normally less than those associated with shortages.
The production quantities are subject to the following capacity restrictions
(6) X'W** 5 ^ for all;,,
and
(7) S 7 /*^ * for all;,,,
where t) is the time required to produce one unit of production stage ;'s product i and Rjt (Ojt), the
capacity in time units of regular time (overtime) available on stage; during period t.
Workforce Considerations
The production decisions can only be carried out if the necessary workforce is available. The
workforce decisions are to be made in view of the production requirements and with respect to the
dynamic costs involved. (For a discussion of dynamic costs see [5]). The workforce available during
period t will be denoted by W t . According to [6] we can differentiate between "productive" (Wf)
and "idle" workforce {W}). So
This distinction is desirable in view of the possibility of keeping workers on a reduced pay rate rathei
than laying them off. Let y) be the number of units of stage ;'s product i which can be produced bj
one worker during period t. Assuming a constant worker productivity (i.e., a constant yj), the require<
productive workforce for regular time or overtime production in t can be expressed by
(9) ^is^ii
196 C. H. von LANZENAUER
and
<io) ri22$*ir
The productive workforce in period t is therefore
(id pf=£i^c^+*jw
The total workforce available in period £ can then be expressed by
(12) £i^(^+^) + ^*i£^(^«+^io) + »7i + ft»^i for all*,
with G?_i being the number of workers hired at the end of period t— 1 and Ff_i the number of workers
laid off at the end of t— 1. Since the decision of which workers to lay off is usually made on grounds of
seniority, we don't have to specifically account for this fact within the model. The same holds for the
"idle" workforce which is generally determined by juniority. The number of workers to be assigned to
the various production stages during period t can easily be obtained from (9) and (10).
From an organizational point of view, it is necessary to maintain the proper relation between the
number of processing workers and supervisory labor. Let Wf be the number of supervisory labor in
period t, then we have an analogous situation to (12)
(13) W t s =W t s . 1 + M l .iN t  1 for all f,
with M t \ being the number of supervisory labor added and N t \ the number laid off at the end of
period t— 1. The relation between the processing and supervisory labor can then be expressed as
(i4) \t±j j V!m+xj»)'*rf.
j=l i
with 8 being the span of control.
The Objective Function
The objective of the model is to determine the production and employment levels so that the con
tribution to profits over the planning period is maximized.
Let r'j t be the selling price of one unit of production stager's product i during t. If we assume that r'
is constant for all t , the total revenues of the production system are
iiirjiDjsj,),
t = i j = i i=i
since any amount backlogged into a future period will bring the same revenue. Assuming an over time
changing selling price rj t , the total revenue depends on the actual price any backlogged demand can
ultimately be sold for. If backlogged demand Bj t has to be filled at the existing price rk of period t the
total revenues are
t=i j=i i=i
MULTISTAGE PRODUCTION SYSTEMS 197
If on the other hand backlogged demand will be filled at current prices, the total revenues become
t=i }=i i=i
Production costs are the sum of material and labor costs. Capital cost will be considered as fixed. The
cost to produce stager's product i depends not only on the decisions made on stagey but also on the
production and employment decisions made on the lowernumbered stages which provide input to
stagey. However, if one handles the labor costs en bloc, the material costs of stagey 's product i are only
the costs for additional raw material since the costs of the inputs from lower numbered stages are
already accounted for. Let cj t be the cost of additional raw material to produce one unit of stage y's
product i during period t. The total material costs are then
t=\ j=i 1=1
The inventory costs are the product of the average amount in inventory and an inventory cost factor
z). Therefore, we have
j
t=i j=i i=i ^
Backlog and shortage costs are
t=\ j=i i=i
with bf being the penalty for backlogging one demand unit of stage /s product i for one period, s! is
is the expected cost in terms of loss of goodwill, over and above the lost revenue, of one unit of stage
/s product i that cannot be filled.
The costs associated with the laborforce are
t =i j=i ,=i J
with p being the regulartime pay rate of a processing worker per period, k the overtime premium,
and A the pay reduction for the idle workforce. h(hs) is the hiring cost per worker (supervisor), f(fs)
the cost of laying off one worker (supervisor), and ps the supervisor's pay rate per period. Combining
the above expressions* leads to the objective function (15)
/
I
t=\ j=\ i=\
(is) tp=z j 2 Wsp^«+*y
"Assuming an overtime constant price r' .
198 C. H. von LANZENAUER
 hG, ~/F t  p s Wf  h s Gf fsFf}
iii{«(4+^(4+*^.
z=i j=\ i=i "■ N ry \ Yj '
 z\{H\ +H i l )b i .B i  (r?+sf)S*
p{\\)WlhGtfF t p s W s t h s Mtf s Nt}
Maximizing (15) with respect to (3) or (3a), (4), (5), (6), (7), (12), (13), (14), and the nonnegativity constraints
for the variables is a problem that can be solved by linear programming.t
III. CONCLUSION
The purpose of model was to formulate the production and employment scheduling problem
on a disaggregate basis in multiproduct, multistage production systems. The model determines the
amount of the demand for each product that should be satisfied, be backlogged, or remain unfilled.
These results are important whenever the production capacity is insufficient to produce all market
requirements. On the other hand, the formulation can be very useful under certain dynamic market
conditions. The model determines also how the productive workforce is to be assigned to the various
production stages. The model seems to be applicable to integrated iron and steel manufacturing systems
as well as to various chemical processes.
The author wishes to thank Andrew Grindlay and Howard E. Thompson for valuable comments and
suggestions.
REFERENCES
f 11 Buffa, E. S., Operations Management — Problems and Models (John Wiley & Sons, Inc., New York, 1968)
[2] Elmaghraby S. E. and A. S. Ginsberg, "A Dynamic Model for the Optimal Loading of Linear MultiOperations Shops,"
Management Technology, Vol. 4, No. 1 (1964).
[3] Haehling von Lanzenauer, C, "Formulation of the Production Scheduling and Sequencing Problem by Bivalent Linear
Programming," Management Science, forthcoming.
[4] Hanssmann, F. and S. W. Hess, A Linear Programming Approach to Production and Employment Scheduling, Management
Technology, Vol., No. 1 (1960).
[5] Holt, C, F. Modigliani, J. Muth, and H. A. Simon, Planning Production, Inventories, and Work Force (PrenticeHall, Englewood
Cliffs, N.J., 1960).
[6] Lippman, S. A., A. J. Rolfe, H. M. Wagner, and J. S. C. Yuan, "Optimal Production Scheduling and Employment Smoothing
with Deterministic Demands," Management Science, Vol. 14, No. 3 (1967).
[7] Silver, E. A., "A Tutorial on Production Smoothing and Workface Balancing," Operations Research, Vol. 15, No. 6 (1967).
[8] Zangwill, W. I., "A Deterministic MultiproductMultifacility Production and Inventory Model," Operations Research, Vol. 14,
No. 3 (1966).
tNote that ^ X 2 ri P\t is a con stant.
A MUTUAL PRIMALDUAL LINEAR PROGRAMMING ALGORITHM
Milton Y. Harris
Naval Research Laboratory
ABSTRACT
A new primaldual linear programming algorithm is exhibited. A proof is given that
optimal solutions to both primal and dual problems (when such solutions exist) are found in
a finite number of steps by this algorithm. A numerical example is included to illustrate
the method.
1.0 INTRODUCTION
Since the publication of the primaldual algorithm by Dantzig, Ford, and Fulkerson, in 1956 (see
Refs. [3], [5], [7]), apparently little work has been done in this area (see Gomory and Balinski [1]). This
paper presents a primaldual algorithm which is similar to that of Dantzig, Ford, and Fulkerson in that
it uses the complementary slackness concept and alternates between primal and dual problems; how
ever, it differs in that it is simpler, but requires a separate Phase I stage.*
Essentially, this algorithm consists in generating and solving a sequence of subproblems which
use the constraints of the original problem. The objective function at each stage, however, is determined
by the solution of the problem of the previous stage. Whether this algorithm is any more efficient than,
say, the revised simplex algorithm remains to be seen.t
Before exhibiting the algorithm, we first establish some notational conventions and state some
necessary results.
2.0 NOTATION AND BASIC THEOREMS
For the sake of simplicity, and without any loss of essential generality, we will consider only
problems of the form
maximize z—c'x
subject to
Ax^b
(2.1)
x^O
where A is an {nxm) matrix, b is an {nx\) vector with nonnegative entries, c is an (mxl) vector of
criterion function coefficients, and x is an {mxl) vector of variables. The "primes" attached to symbols
denote the transpose of the indicated matrix. All entries in /I, b, and c are taken to be given constants.
We may state the dual to this problem as
(2.2) minimize u—w' b
*For an explanation of the two stages of the simplex method, see Dantzig [4], ch. 5.
t A Fortran program has been written for the CDC 3800 computer to compare the efficiency of this algorithm with the sim
plex method, but as yet, no conclusive results have been obtained.
199
200 M. Y. HARRIS
subject to
w'A S= c'
(2.3)
w' ^0
where w is an (nxl) vector of variables.
These two problems may be rewritten by introducing slack variables 5 and v for the primal and
dual, respectively, and taking the transpose of 2.3 as follows:
maximize z — c'x minimize u — w' b
subject to subject to
Ax + Is = b A'w — Iv = c
(2.4) (2.5)
{x,s)^0 {w,v)^0,
where / is an appropriately dimensioned identity matrix in each case.
We will digress here slightly to define formally the notion of complementary slackness, a concept
which is central to the algorithm.
DEFINITION: For given basic feasible solutions (x, s) to the primal and (w, v) to the dual,
complementary slackness is satisfied if, and only if,
i) w • s —
and
ii) x  v — 0.
We conclude this section by stating four wellknown results which are essential to our algorithm.
Proofs are given by Dantzig in [4].
Suppose we have any two basic feasible solutions (x, s) for the primal (2.4) and (w, v) for the
dual (2.5). Then
LEMMA 1: w'b^c'x.
LEMMA 2: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and
only if, w'b = c'x.
THEOREM 1: (x, s) and (w, v) are optimal solutions for the primal and dual, respectively, if, and
only if," complementary slackness is satisfied for these solutions.
FUNDAMENTAL DUALITY THEOREM: If a feasible solution to 2.4 exists and z of 2.4 is bounded
above, then optimal feasible solutions exist for both the primal (2.4) and the dual (2.5).
3.0 THE ALGORITHM
The algorithm is initiated with basic feasible solutions to both the primal and dual, either or both
of which may be contained in the intital formulation of the problem or may be generated by Phase I of
the simplex algorithm. We also assume that z of 2.4 has a finite upper bound. The fundamental duality
theorem assures us that optimal solutions for both the primal and dual exist.
Recall the primal and dual constraints may be written as
(3.1) Ax + Is = b (primal)
and
(3.2) w'Av'I=c' (dual).
MUTUAL PRIMALDUAL ALGORITHM . 201
Multiplying 3.1 by w' on the left and 3.2 by x on the right and applying Lemma 1 gives
(3 3) ws + vx = wb — cx^0.
Let a = w • s + v • x. Now any decrease in the value of a brings us closer to the optimum since a decrease
in a represents a movement toward the satisfaction of complementary slackness and a convergence of
wb and cx toward their common optimal value (see Lemma 2 and Theorem 1).
We now proceed to define the algorithm. Suppose we have a basic feasible solution {x°, s°) to the
primal and a basic feasible solution (w°, v°) to the dual with
(x°, 5°) • {v°, w°)=x v + s°w >
(i.e., the solutions are nonoptimal). Let Y° be a current fulltableau matrix for the primal, F° be a cur
rent constant column for the primal, Z 1 be a current fulltableau matrix for the dual, and Z* be a current
constant column for the dual. Consider the linear programming problem
(3.4) minimize a — (x, s) • (v°, w°)
subject to
(3.5) Y°(x,s) = Y°
(x, s) ^ 0.
Suppose the optimal (if one exists) solution to 3.4 and 3.5 is (x 1 , s 1 ) and
(zSs 1 )^ ,^ )^.
Let Y 1 and YJ be the matrix and constant column, respectively, corresponding to this solution. If
6*^ = 0, then {x 1 , s 1 ) and (w°, v°) are optimal solutions to the primal and dual, respectively, by Theorem
1. If oto > (recall from 3.3 that a<f 3= 0), consider the problem
(3.6) minimize a^ (x l , s 1 ) • (v, w)
subject to
(3.7) Z*{w i v)=Z\
(W, V) 52=0.
Let the optimal tableau to this problem be represented by Z 2 and Z\ and the optimal solution be
(w 1 , v 1 ); let
(x^s 1 ) ■ (v 1 ,w 1 )=at
As before, if af = 0, we are through. If not, we return to the primal and solve
(3.8) minimize a 2 = (x, s) • (v 1 , w l )
subject to
(3.9) °
(x,s) >0
for (x 2 , s 2 ), Y 2 , Y§, and a* = (x 2 , s 2 ) ■ (v 1 , w x ). The algorithm consists in following this procedure of
using the optimal solution at each stage as the criterion function coefficients for the next stage until
for some n, a* = 0.
202 M. Y. HARRIS
In the next section we will show that the algorithm results in optimal solutions for both the primal
and dual, whenever they exist, in a finite number of iterations.
4.0 PROOF OF CONVERGENCE
We may state the subproblems described above in general as follows:
Primal
minimize ak = (x, s) ■ (v^ k \ w [k] )
subject to
YM(x,s) = YM
(4.1)
(*, s)2=0,
where [k] — k/2 for A; = 0, 2, 4, . . . and
Dual
minimize a m = (x lm \ s [m] ) ■ {v, w)
subject to
ZM(w, v)=ZW
(4.2)
(w, v) 5=0,
where [m] = (m— l)/2+ 1 for m — 1, 3, 5, . . . We first show that the subproblems 4.1 and 4.2 have
optimal solutions.
THEOREM 2: Both problems (4.1 and 4.2 above) have optimal solutions for any ke{0, 2, 4, . . .}
and any me{l, 3, 5, . . .}.
PROOF: Both problems (4.1 and 4.2) have feasible solutions by assumption. Now for every basic
feasible solution (x, s) to the primal and (w, v) to the dual, we have by Lemma 1
(4.3) x v + w s^O.
Since the problems in question are both minimizing problems whose objective functions are of the
form of 4.3, we have just seen that both are bounded below by zero. By the fundamental duality theorem
then, both have optimal solutions.
Finally, we show that the algorithm does indeed converge.
THEOREM 3: a„* = for some (finite) n.
PROOF: We assume that the original problem has an optimal solution. Clearly
a* < ao,
since at worst (w°, v°) satisfies 3.7 and (x 1 , s 1 ) • (v°, w°) — c*o (this fact is just a general property of
constrained optima, that is, that the optimal solution to a minimizing problem is no larger than any
admissable solution). Suppose a* = a*. Then for every basic feasible solution (w, v) for the dual
w • b — c • x x = w • s x + s l + x l • v (by 3.3)
= «i? a* = ctQ (assumption)
= w° • s 1 + x 1 ■ v°
= w°bcx 1 (by 3.3).
MUTUAL PRIMALDUAL ALGORITHM 203
Thus for every basic feasible solution (w, v) for the dual
w • b — c • x 1 2= w° • b — c • x 1 or w • b 2= w° • b,
which implies that (w°, v°) is the optimal solution for the dual. Similarly, for every basic feasible
solution (x, s) for the primal
w° • b — c • x = w° • s + x • v° (by 3.3)
= «i 2= a* — ct * (assumption)
= w°bcx 1 (by 3.3)
Thus for every basic feasible solution (x, s) for the primal
w° • b — c • x ^ w° • b — c ■ x l or c • x ^ c ■ x 1
which implies that (x l , s 1 ) is optimal for the primal. Now if a* > 0, then complementary slackness
is not satisfied by the solutions (*', s 1 ) and (w°, v°) , and so at least one must be nonoptimal by Theorem
1. This contradiction (i.e., we just showed that if a* — a*, both (m; , v°) and (x 1 , s 1 ) are optimal) im
plies that if a* > 0, then
a* < a *.
We may repeat this argument for ai and ocz and so on. We conclude therefore that
a * > implies that a* +1 < a* for each A; = 0, 1,2,. . ..
Since there are only a finite number of basic feasible solutions to the primal and dual, and since
each iteration of the algorithm results in a strict decrease in a, if optimal solutions exist, then by
Theorem 1, there exists an n (finite) such that
Furthermore the solution to the problem for which a n is the objective function is the optimal solution to
either the primal or dual with the solution to the preceeding problem being the optimal solution to the
other.
5.0 NUMERICAL EXAMPLE
Consider the following primal problem (this problem was taken from Hadley [5], p. 269):
maximize z= 2x x + 13* 2 + Tx 3 + 6x4
subject to
1 & 4*i + lx 2 + 3*3 + 5x 4
3 3= — xi + 3x2 + 3x 3 — 3x 4
(5.1) 7 3*2*1 4x 2 +lx 3 + 8x4
5 5s 3xi + 5x 2 — 2x 3 — 4x 4
2 2= — 5xi — lx 2 + lx 3 + 6X4
8 5=6xi + 2x29x3+lx4
The dual is then
minimize u = W\ + 3k;2 + 7w 3 + 5w 4 + 2w 5 + 8w 6
204
M. Y. HARRIS
subject to
(5.2)
2 =£ 4wi — 1^2 + 2tt>3 + 3u>4 — 5w$ 4 6wa
13 =S 1 w\ + 3z#2 — 4m; 3 4 5t^4 — \w$ + 2we
7 =S 3wi + 7w 2 + Iwa — 2w4 + lws — 9w 6
6 =£ 5^i — 3^3 + 8w 3 — 4w 4 + 6^5 + 1m;6
After appending slack variables S\ through se to 5.1 and v\ through y 4 to 5.2 and writing both systems
in detached coefficient form, we have
Primal
(5.3)
b
Xi
*2
x 3
*4
Sl
«2
53
54
55
56
2
13
7
6
1
4
1
3
5
1
3
1
3
7
3
1
7
2
4
1
8
1
5
3
5
2
4
1
2
5
1
1
6
1
8
6
2
9
1
1
(5.4)
Dual
c
W\
Wi
w 3
Wi
w 5
w$
fli
V2
v 3
Va
1
3
7
5
2
8
2
4
1
2
3
5
6
1
13
1
3
4
5
1
2
1
7
3
7
1
2
1
9
1
6
5
3
8
4
6
1
1
Carrying out a phase / step of the simplex algorithm on 5.4 gives us the dual tableau (see Table 5.5)
in basic feasible form (the objective row has been omitted since it is not needed in the rest of the
procedure).
Now we construct an objective row for 5.3 from the basic feasible solution to the dual (5.5). This
is shown in Table 5.6.
After adjoining 5.6 to 5.3, four pivots of the simplex algorithm are required to minimize. The result
is shown in Table 5.7.
Since the objective function has not been reduced to zero (entry (1, 1) in 5.7), at least one dual
cycle is required. We construct an objective row for the dual (5.5) from the solution to the primal (5.7).
This is shown in Table 5.8.
After adjoining 5.8 to 5.5, one pivot of the simplex algorithm is required to minimize. The result
is given in Table 5.9.
Since the value of the objective function in 5.9 is now zero, 5.9 must be the optimal dual tableau
(except for the objective row, of course). Thus the previous primal tableau (5.7) must have been optimal.
Checking the values of the actual objective functions using the solutions given in 5.7 and 5.9 for the
primal and dual, respectively, we find that both are 13.
MUTUAL PRIMALDUAL /LGORITHM
205
ACKNOWLEDGEMENTS
The author wishes to thank Dr. Richard Young, Associate Professor of Economics, Rice University,
and Mr. Arthur Ziffer of the U.S. Naval Research Laboratory for their many helpful suggestions toward
the development of this algorithm.
Table 5.5
c
w,
W2
w 3
w 4
W 5
We
V\
V2
Va
v 4
1.2979
1.0000
0.0000
0.8511
0.0000
0.0000
0.3830
0.0164
0.0213
0.0426
0.0851
2.2098
0.0000
0.0000
0.7328
1.0000
0.0000
1.0040
0.0247
0.2044
0.0797
0.0267
0.8179
0.0000
1.0000
0.4206
0.0000
0.0000
1.1544
0.0351
0.0302
0.1069
0.0421
1.8006
0.0000
0.0000
0.0747
0.0000
1.0000
0.0604
0.1227
0.1336
0.0351
0.0925
Table 5.6
b
Xl
x 2
X3
Xi
5l
52
5 3
Si
5 5
5 6
18.4018
0.0000
0.0000
0.0000
0.0000
1.2979
0.8179
0.0000
2.2098
1.8006
0.0000
Table 5.7
6
Xi
X 2
x 3
X4
5l
S2
S3
S4
$5
56
5.4018
0.0000
0.0000
51.2537
0.0000
3.4789
6.2958
0.0000
1.8733
0.0000
0.0000
0.0000
1.0000
0.0000
3.4930
0.0000
0.0423
0.4085
0.0000
0.2535
0.0000
0.0000
.0.0000
0.0000
0.0000
2.6338
1.0000
0.1972
0.2394
0.0000
0.1831
0.0000
0.0000
11.0000
0.0000
0.0000
2.1268
0.0000
0.7606
0.6479
1.0000
0.5634
0.0000
0.0000
1.0000
0.0000
1.0000
3.8028
0.0000
0.1831
0.4366
0.0000
0.0986
0.0000
0.0000
3.0000
0.0000
0.0000
28.4648
0.0000
1.2113
3.0423
0.0000
2.2676
1.0000
0.0000
6.0000
0.0000
0.0000
1.7183
0.0000
0.3099
1.3380
0.0000
1.1408
0.0000
1.0000
Table 5.8
c
Wi
U>2
W 3
Wa
M>5
We
V,
V2
va
V4
5.4018
0.0000
0.0000
11.0000
0.0000
3.0000
6.0000
0.0000
1.0000
0.0000
0.0000
Table 5.9
c
W\
Ul2
w 3
W4
(#5
We
Vl
V2
Va
V4
0.0000
0.0000
0.0000
11.0000
0.0000
3.0000
6.0000
0.0000
1.0000
0.0000
0.0000
2.8593
1.0000
0.0000
0.7863
0.0000
0.8672
0.3306
0.0000
0.1372
0.0122
0.1653
1.8473
0.0000
0.0000
0.7178
1.0000
0.2013
1.0162
0.0000
0.1775
0.0726
0.0081
0.3028
0.0000
1.0000
0.3992
0.0000
0.2861
1.1371
0.0000
0.0080
0.1169
0.0686
14.6748
0.0000
0.0000
0.6088
0.0000
8.1500
0.4923
1.0000
1.0888
0.2861
0.7539
REFERENCES
[1] Balinski. M. L. and R. E. Gomory, "A Mutual PrimalDual Simplex Method." in Recent Advances in Mathematical Program
ming (R. L. Graves and P. Wolfe, editors) (McGrawHill Book Co., New York, 1963), pp. 1726.
[2] Dantzig, G. B., "Composite SimplexDual Simplex AlgorithmI," The RAND Corp., RM1274 (Apr; 1954).
[3] Dantzig, G. B., L. Ford, and D. Fulkerson, "A Primal Dual Algorithm for Linear Programs," Paper 7 in "Linear Inequalities
and Related Systems," Annals of Mathematical Studies (Harold Kuhn and Albert Tucker, editors) (Princeton University
Press, Princeton, 1956), No. 38, pp. 171181.
[4] Dantzig, G. B., Linear Programming and Extensions (Princeton University Press, Princeton, 1963), Ch. 5 and 6.
206 M  Y  HARRIS
[5] Hadley, G., Linear Programming (Addison Wesley, Reading, Mass., 1962), pp. 257266.
[6] OrchardHays, W., "A CompositeSimplex AlgorithmII," The RAND Corp., RM1275 (May 1954).
[7] Vajda, S., Mathematical Programming (AddisonWesley, Reading, Mass., 1961), pp. 108113.
INTERCEPTION IN A NETWORK*
R. D. Wollmer
The RAND Corporation
Santa Monica, California
ABSTRACT
This paper presents an algorithm for determining where to place intercepting units in
order to maximize the probability of preventing an opposing force from proceeding from one
particular node in an undirected network to another.
The usual gaming assumptions are invoked; namely, the strategy for placing the units
is known to the opponent and he will choose a path through the network which, based on
this knowledge, maximizes his probability of successful traverse. As given quantities, the
model requires a list of the arcs and nodes of the network, the number of intercepting units
available to stop the opposing force, and the probabilities for stopping the opposition at the
arcs and nodes as functions of the number of intercepting units placed there. From these
quantities, the algorithm calculates the probabilities for placing the unit at the arcs and nodes
when one intercepting unit is available, and the expected numbers of units to place at the
arcs and nodes when multiple intercepting units are available.
I. INTRODUCTIOIV
The paper depicts a game between two opposing forces — an evader and an interceptor. Only the
interceptor's problem will be dealt with. In this game, the evader attempts to proceed from one point
to another in a network. The interceptor, who may possess one or more indivisible intercepting units,
attempts to stop him by placing these units along arcs or nodes that he expects the evader to traverse.
His problem is to place these so as to minimize the evader's probability of successful traverse. The
evader's problem is, of course, to select a path that will maximize this probability. Both problems can
be represented by a zerosum twoperson game. As the next section will show, however, the game matrix
is very large and difficult to generate. By using a steepest ascent approach, a solution may be obtained
much more easily for the interceptor.
This paper presents an algorithm based on such an approach. As inputs, it requires a list of the
arcs and nodes of the network, the number of intercepting units available to stop the opposing force,
and the probabilities for stopping the opposition at the arcs and nodes as functions of the number
(integer) of intercepting units placed there. From this information, the algorithm calculates the expected
number of intercepting units to place at each arc and each node. For the case where only one inter
cepting unit is available, the expected values are probabilities of placement; for more than one intercept
ing unit, expectations are to be interpreted in the most obvious way. Specifically, if the expected number
of units to place at a particular location is 41/4, one would always place four there and would place
a fifth one there onefourth of the time.
The situation depicted by this algorithm is particularly relevant in applying interdiction to infiltra
tion and counterinsurgency, where the counterinsurgent or infiltrating force would play the role of the
*This research is supported by the United States Air Force under Project RAND. Views or conclusions contained in this
study should not be interpreted as representing the official opinion or policy of the United States Air Force.
207
208 R D. WOLLMER
evader, while the interdiction forces would play that of the interceptor. It also applies to police networks
where a wanted fugitive (evader) is attempting to escape from a city while police (intercepting) units
are attempting to close off his avenues of escape.
The algorithm to be presented is a steepest ascent approach and yields solutions that are always
optimal in the gaming sense (i.e., the evader's best chance of successful traverse is as small as possible)
for the interceptor when he has but one intercepting unit at his disposal. For the case of multiple inter
cepting units, the solutions are not always optimal; thus for the general case the algorithm is an
approximation scheme.
II. INTERCEPTOR'S PROBLEM FOR ONE INTERCEPTING UNIT
The network is characterized by sets of elements called arcs and nodes. Nodes are points or
junctions, and arcs are undirected line segments joining nodes. An arc that joins node i to node j may
be disignated by the symbol (i,j) or (J, i). The node from which the evader starts is called the source;
the one he attempts to reach is called the sink.
It is assumed that the evader and interceptor are each composed of one indivisible unit. Later on
this restriction will be relaxed for the interceptor. The evader is attempting to choose a path from
source to sink in such a way that his probability of successful traverse is maximized. The interceptor
wishes to place his intercepting unit so that this probability is minimized. As is customary in game
theory, it is assumed that each side's strategy will become known to the other and consequently each
will need to adopt a mixed strategy. The inadequacy of pure strategies for the interceptor may be dem
onstrated by any network whose source and sink are relatively invulnerable to attack and any particular
intermediate arc or node may be bypassed. Placing the intercepting unit at the source or sink would be
relatively ineffective while placing it anywhere else would be completely ineffective as the evader
would choose a path that bypassed the unit.
Quantities defining the local effectiveness of placing the intercepting unit at arcs and nodes,
and quantities expressing the interceptor's strategy are as follows:
p(i) = probability that the evader will be stopped at node i, given that node i is on his chosen path
and the interceptor places his intercepting unit there.
p(i, j) = probability that the evader will be stopped at arc (i, j) given that arc (i, j) is on his
chosen path and the interceptor places his intercepting unit there.
77 (i) — probability that the interceptor places his intercepting unit at node i.
7r(i,y) = probability that the interceptor places his intercepting unit at arc (i,j).
The probability that the evader can successfully cross a particular node, t, is 1 — n(i)p(i) and
that of crossing a particular arc (i,j) is 1 — Tr(i,j)p{iJ). Assuming these probabilities are independent*
and the evader attempts to reach the sink by the path i,, . . . i n , his probability of successful traverse is
n n1
(1) ^ = E[ [ l1T (ir)p(ir)] ][ [1— ir(ir, ir+l)p(ir, U+l)]
r=l r=\
This quantity, K, will be referred to as the value of the path ii, . . ., l».
Given the p{i), p(i,j), tt(i), and v(i,j), the evader's problem is to find a sourcesink path of
maximum value. The interceptor's problem, given the p(i) and p(i, j), is to choose 7r(t) 2*0 and
*This assumption is in fact not true. To illustrate this, suppose ir(i) = l/2 and n(j) = 1/2. Independence would require
that the probability of the intercepting unit being at nodes i and j at the same time be 1/4 while in reality it is zero. However,
the next section will show that this assumption leads to no inaccuracies for the oneinterceptingunit case.
INTERCEPTION IN A NETWORK ._ 209
7r(i, j) 3 s such that the maximum value of all sourcesink paths is minimized, subject to the con
straint %Tr(i) +2,ir(i,j) *£ 1. Note that if all p(i) and p(i,j) are strictly positive, any optimal strategy
will be such that strict equality holds in this relationship. To see this, suppose the intercepting unit
is not placed at any arc or node with probability p (i.e., 2/7r(t) + 2/7r(i, j) = 1 — p). Then increasing
all ir(i) and Tr(i,j) by equal amounts totaling p yields an improved solution.
III. A STEEPEST ASCENT APPROACH FOR ONE INTERCEPTING UNIT
The problem of finding the n{i) and ir(i,j) can in theory be solved by linear programming. The
linear program solution would include an optimal strategy for the evader as well as the interceptor.
The value of the game would be equal to both the evader's maximum guaranteed probability of
successful traverse, and one minus the interceptor's guaranteed probability of stopping the infiltrator.
However, the game matrix for the program would require a column for every possible sourcesink
path and a row for each of the interceptor's pure strategies. The number of sourcesink paths is exceed
ingly large even for relatively small networks and in later sections, when multiple intercepting units
are allowed, the same may be said of the number of pure interceptor strategies.
The algorithm of this section is a steepest ascent approach which avoids these enumeration
problems; however, while it yields an optimal strategy for the interceptor and a value for the game, it
does not yield an optimal strategy for the evader.
Initially, all 7r(i) and ir{i,j) are assigned values of zero. Then they are increased by small amounts
A7r(i) 2=0 and A7r(i, j) 2=0. The proportions for the Air(i) and ^rr(i,j) are such that the additional
intercepting unit placement probability, 2A7r(£) +XA7r(i, _/), divided by the decrease in maximum
path value, is minimized as these two quantities tend toward zero. This is equivalent to maximizing
the decrease in maximum path value per unit of additional intercepting unit probability allocation.
The 7r(t) and ir{i, j) are then increased again in the same manner until they sum up to one.
For a specific path, t'i, . . ., i n , increasing Tr(i r ) by A7r(i r ) replaces the factor [1 — ir(i r )p(i r )]
in Eq. (1) by [1 — (rr(i r ) + A7r(£ r ))p(£ r )J and reduces the value of the path from K to K — &K, where
1 lir(»r)p(»r) J"
Solving for the quotient of the additional intercepting unit probability allocation and the decrease in
path value, the following expression is obtained:
(6a) AK K
ir(ir)
P(lr)
If intercepting unit placement probability were increased at an arc instead of a node, the same
argument would yield an expression for the additional intercepting unit probability placement divided
by the decrease in path value which is identical to the preceding one, with arc probabilities substituted
for node probabilities. Specifically, for arc (i r , ir+i), this would be
(3b) **&' ir+l) ~
7T(i r , Ir+l)
A£ Klp(i r ,i r+1 )
Note that in order to obtain any decrease in maximum path value, the A7r(£) and kir{i,j) must be
strictly greater than zero along a subset, C, of arcs and nodes that interesects all paths of maximum
210 R D. WOLLMER
value. An efficient allocation would require that all members of C remain on maximum value paths,
for otherwise, lessening the probability increase at an arc or node not on a maximum value path by a
small amount would decrease SA7r(i) + 2A7r(i, j) , while not affecting the maximum path value. Then,
of course, this small amount of probability could be redistributed among arcs and nodes on maximum
value paths to get a strict improvement. Thus, the Air(j) and krr(i,j) assigned to the chosen C must be
such that the A£'s are equal for each of its members. In other words, they must be proportional to
77(0
(4) p(i)
b) , 7T(i,j).
p(l,J)
The constant — , where K is the maximum path value, was dropped in Expression (4), of course. The
total increase in intercepting unit placement probability per unit decrease in the value of a maximum
probability path would be:
n(i,j)
p(ij)
Thus, the problem is reduced to finding a subset of arcs and nodes intersecting all maximum prob
ability paths that minimizes Expression (5).
Initially, all paths have value one since all ir(i) and 7r{i, j) are zero. Thus C is required to be a
set blocking all sourcesink paths. If the nodes and arcs are assigned capacities equal to the quantities
of Expression (4a) and (4b), the problem of finding C is reduced to one of finding a minimum cut.*
However, this is equal to the value of the maximum sourcesink flow and can be found by the maximum
flow algorithm of Ref. [3].
Note that once additional intercepting unit probability is allocated to the minimum cut, C, the
Capacity of each arc and node in C decreases by its additional probability increase while the capacity
of any other arc or node remains the same. Thus, the value of C will be reduced by SA7r(i) +SA7r(i, j)
while all other cuts will decrease by amounts that do not exceed this. Hence, C remains minimum and
the entire unit of probability may be allocated to C.
The 7r(i) and ir{i, j) may therefore be solved for as follows:
1. Assign all nodes capacities of llp(i) and all arcs capacities of l/p(i, j) (since all 7r(i) and
7r(i, j) are zero).
2. Maximize flow from source to sink.
3. Let C be the minimum cut set and ^ its value. Set ir(i) =p(i)/V for ieC and Tr(i,j) =p(i,j)IV
for (i,j)eC. Set all other 7r(i) and Tr(i,j) equal to zero.
Note that since only one cut is obtained, each path of maximum value has positive intercepting
unit probability allocated to only one of its arcs or nodes, and both its value and traverse probability
are equal to the probability associated with that particular arc or node. Hence, the independence
assumption does not lead to inaccuracies for the oneinterceptingunit case. Finally, the maximum path
value is 1 — p, where V is the value of C found in step 3.
*Essentially, a cut set is a set of arcs and nodes blocking all sourcesink paths. Its value is the sum of the capacities of its
arcs and nodes, and the minimum value of all cuts is equal to the maximum value of the sourcesink flow. See Ref. [3 J.
INTERCEPTION IN A NETWORK 211
IV. MULTIPLE INTERCEPTING UNITS
The steepest ascent algorithm of the last section may be modified to handle the case where there
is more than one intercepting unit; however, the interceptor strategy obtained in some cases may not
be optimal and the game value obtained may be approximate rather than exact. This will be discussed
more fully later on. The quantities measuring the local effectiveness of placing intercepting units at
arcs and nodes and those defining the interceptor's strategy must be redefined as follows:
p(i)k — probability that the evader will be stopped at node i, given that node i is on his chosen
path and the interceptor places k intercepting units there.
pi.ii J )k — probability that the evader will be stopped at arc (i,j), given that arc (i,j) is on his chosen
path and the interceptor places k intercepting units there.
tt{i) = expected number of intercepting units the interceptor places at node i.
ir(i,j) — expected number of intercepting units the interceptor places at arc (i,j),
where p(i)o andp(i,j)o are both identically zero. Also let
7r(i')=7r(;) — [«■(*)]
ir{i,j) =ir(i, ;') — [n(i, j) ] .
Note that 7r(i) and 7r(i,j) are merely the fractional parts of 7r(i) and n(i,j). It will be assumed
p(i)k and p(i,j)k are increasing in k and
(6) a) p(i)k+i—p(i)k^p{i)k—p(i)ki
b) p(i, j)k+i—p(i, j)k < p(i, j)k~ p{i, j)ku k > 1.
The assumption of Expression (6) essentially states that if an additional intercepting unit is placed
at any arc or node, it cannot be more effective in stopping the evader than any previous unit placed
there. It is often satisfied in practice as the effects of overkill become more pronounced as the number
of units at a location increase. The simplest example of this is the case where each unit placed at node i
has independent probability p of stopping the evader. For this situation p(i)fc— 1 ~~ (1 — p) k and
p(i)k+i~ p(i)k = p(l— p) k which is certainly decreasing in k, thus satisfy (6).
The values of 7r(i) and n(i,j) of course are not sufficient in themselves to define the interceptor's
strategy as any nonzero value of 7r(i) or ir(i,j) may often be realized in several ways. However, the
assumption of Expression (6) assures that for a given value of n ( i ) , the probability of stopping the evader
at node i should he attempt to cross it is minimized by allocating [7r(t)] + 1 units at i with probability
7r(i), and [7r(i)] units with probability 1 — Tr(i). In other words, if 7r(i) = 31/2, tt{i) would best be
realized by assigning three forces at i half the time and four half the time, as opposed to such a policy
as allocating two forces half the time and five half the time. This is shown in the theorem to follow. First,
define P{k, i} as the probability that k units are placed at node i. Note that 2,kP{k, i} = ir(i) and the
probability of thwarting an evader's attempt to cross node i is 1 — 2,P{k, i}p(i)k. Minimizing this prob
ability is equivalent to maximizing 2,P{k, i}p(i)k.
THEOREM: The quantity (t')2P{/:, i}p(i) k is maximized, subject to the constraint (ii)^,kP{k, i}
= ir(i) when P{[ir(i)] + 1, i}=n(0 ,<P{[n(i)], i} = ln(i), and all other P{k, i} = 0.
PROOF: Let h = max {k/P{k, i}>0}, /= min {k/P{k, i} > 0}, and d=hl. Let P{k, i} = P{k, i}
be a set of values of smallest d which maximizes (i) subject to (ii). Suppose the theorem is false.
Then hl = d^2. Let m = min{P{h, i}, P{1, i}}. Set P{h, i} = P{h, i}  m and P{h1, i}
= P{h — l, i} + m. Then set P{1, i} = P{l, i} — m and increase P{/f 1, i} by m. (If h — l = 2, then
P{l+l} = P{l+l} + 2m.) Either P{h, i} or P{1, i} is now zero, decreasing d by at least one unit.
212 R D. WOLLMER
The increase in (/) is equal to m(— p(i)i + p(i)i+\ — p(i)ft + p(0/ii) Since (ii) is still satisfied, this
is strictly negative from our choice of P{k, i}. However, since h > / + 1,
P(i)h — p(i)hi =S p(i)i+i—p(i)i
 p(i)i + p(i)i+i— p(i)h + p(i)hi 2=
for a contradiction. Q.E.D.
The same results, of course, hold for the Tr(i, j) of the arcs. The significance of this is that the
search for an optimal interceptor strategy may be confined to strategies of the above type. When this
is done, the 77(1) and Tr(i, j) are sufficient to specify the interceptor's strategy.
As in Eq. (1) of Sec. II, if independence is assumed, the value of a path is still the product of the
probabilities associated with its nodes and arcs; however, the probability associated with node i is now
7T(j) (1— p(i)[ir(l)]+l) + (117(0) (1— P(i)[n(l)])
or
^~p(i)[Mi)] — rr(i)(p(i)i7T(i)]+i—p(i)[Tr{i)]),
and that with arc (i, j) is
l^P(^y)[7r(i,j)] — 7r(l, j)(p(i, j)[7r(i,j)] + l— P(i, j)[n(i,j)]).
If ti, . . ., in is a path of value K, and 7r(i r ) is increased by A7r(i r ), reducing the path value to
K — &K, the expressions for K and K — &K, as before, differ only in the factor for node i r . Thus, for
sufficiently small A7r(i r ), the expression for the quotient of the intercepting unit increase and path
decrease is
A7T(l r ) 1 f 1—P{ir)[n(i r )) 1
TT{l r ) \ ,
AX K [p(ir)[n(i r )]+l— P(ir)[7r(i r )]
and decreasing path value by increasing the number of intercepting units at an arc instead of a node
yields the analogous expression,
A7T(ir, ir+l) If 1 — P(ir, ir+l)[n(i r ,i r+l ))
A£
If 1 — PKlr, lr + l)[ir(i r ,i r+l )] ... 1
= p\ — — : : T. — : ^ ?rUr, ir+\) \
A I p{l r , lr+l)ln(i r ,i r + 1 )] + l—PKlr, lr+l) [n(i r , i r + 1 )] J
As in Sec. Ill, intercepting units must be allocated along all arcs of some cut set, and the problem
is reduced to finding a minimum cut set, which in turn reduces to one of finding a maximum source
sink flow. The node and arc capacities are now
(7) a) C(l)=—7^. 7TT 77 (I)
P(l)[77(t)]+1P(l)[7r(i)]
b) C(l,J)= ,. ., /• x TT\l,J)
Note that the c(i) and c{i, j) are no longer decreasing functions of 77(1) and 77(1,7) except over
intervals whose endpoints possess the same integer part. Specifically, c(i) is a decreasing function of
INTERCEPTION IN A NETWORK 213
7r(t) for < n(i) < 1, 1 s£ 7r(i) < 2, etc. However, c(i) is discontinuous and actually increases at all
integer points for which (6a) holds with strict inequality. This follows since
lim l—p(t) [„(!)] . ,., _ 1— ■p(Q[ir(i)]+i < 1— p(t)[7r(»)]+i
ff(i) ^ 1 p(i)[n(i)]+l— P(i)[ir(i)] P(i)[n(i)]+1—P(i)[n(i)] P (0 [n(i)]+2 ~ P (0 [ir(i)]+l
Similar results hold for c{i, j). Thus, when a minimum cut set is found and additional intercepting
units are allocated along its arcs and nodes, that cut remains minimum provided the expected number
of intercepting units on these arcs and nodes do not reach or exceed their next highest integer values.
This requires that the steepest ascent approach be modified as follows.
Starting with all 7r(t) and ir(i, j) equal to zero, capacities are assigned to the arcs and nodes of
the network which are equal to the quantities defined by Expression (7). Flow is maximized from
source to sink and the 7r(t) and Tr(i,j) corresponding to the arcs and nodes of the minimum cut set are
increased by amounts A7r(i) and A7r(i, j) which are proportional to their capacities. The constant of
proportionality, M, is the smallest possible such constant that will either increase some 7r(t) or Tr(i,j)
to the next highest integer or will increase %Tr(i) + ?.TT{i,j) to n. If the latter does not happen, the new
values of ir(i) and rr(i,j) are used to calculate new capacities and flow is maximized again to obtain
a new cut. The process is repeated until %7r(i) +£77 (1,7) = n.
Specifically, the algorithm for the general case of n intercepting units is as follows:
1. Set all 7r(i) = 0and all ir{i,j)=0.
2. Assign the nodes and arcs of the network capacities of c(i) and c(i,j) respectively, where c(i)
and c(i, j) are as defined by Expression (7).
3. Maximize flow from source to sink and let C be the nodes and arcs of the minimum cut set found.
A n Tit {1 _ tHO}
4. Compute M x = mm  1 ...
i€C c(i)
1VI2 — mm — —
(i,j*C C(l,j)
5. Let M=min (M u M 2 , M 3 ). Set bn(i)=Mc(i) if ieC and &ir(i,j)=Mc(iJ) if (i,j)eC. Set all
other A7r(t) and Air(i,j) equal to zero.
6. Set ir(i) = ir{i) + Att({) and ir(ij) = n(i,j) + bir(i,j) for all i and (i,j). If %ir(i) = £w<i,» < n,
go back to step 2. Otherwise terminate.
Step 1 is simply an initialization step. Steps 2 and 3 determine which nodes and arcs should have
their expected number of intercepting units increased and the proportions of the increase. Specifically
the expected number of intercepting units will be increased along the nodes and arcs of the minimum
cut set found in proportion to their capacities. Steps 4 through 6 determine the proportionality constant
M, makes the appropriate increases, and makes a termination test. Note that if M = M\, some rr{i)
will be increased by 1 — Tr(i), and since 77(1) is the fractional part of 7r(i), this will increase 7r(i) to an
integer value. Similarly if M=?M 2 , some v(i, j) will be increased to an integer value. If M = M 3 , total
expected intercepting units will be increased by nXTr(i) —%TT(i,j) to n. In the latter case all n inter
cepting units have been allocated and the procedure terminates.
214
R. D. WOLLMER
At termination the ir{i) and rr{i, j) represent the expected number of intercepting units to place
at node i and arc (i, j), respectively. As mentioned earlier, the strategy will be to place [77(1)] + 1
forces at node i with probability ir(i) and [77(1)] forces with probability 1 — Tr(i).
The maximum path value can be kept track of during the course of the algorithm. Specifically,
if K is the maximum path value at the beginning of an iteration, then the maximum path value at the
end of the iteration is K(l — M), where M is that found in Step 5. Of course, K= 1 at the beginning of
the first iteration. This is easily proven by letting V be the value of C, the minimum cut set. Then
substitute the generalized capacities (Expression 7a— b) for the specialized ones (Expression 4a— b)
in Expressions (3) and (5). These replacements for Expressions (3) and (5) will have the same meaning
in the multiple intercepting unit situation as Expressions (3) and (5) did for the single intercepting
unit situation. The replacement for (5), which represents the total increase in intercepting unit place
ment probability per unit decrease in the value of a maximum probability path, is:
^Air(i) +SAn(t, /) = 1 y
AK K
However, 2A7r(i) +£A7r(i, j) =MV where M is as defined in Step 5 of the multiple intercepting unit
algorithm. Thus
MV = \y
AK K
or
KAK = K(lM).
Note that when there is more than one intercepting unit to be assigned, the strategy obtained is
not necessarily optimal. To illustrate this, consider the network of Fig. 1.
Figure 1. A simple network with node 1 the source and node 6 the sink.
Suppose p(l, 2)i = 5/8, p(l, 2)2=1, p(2) fc = l — (1/4)*, and all other p(i)k and p(i, j) k equal
1 — (1/2)*. If there were two intercepting units available, placing both at arc (1, 2) with probability
one will make the evaders successful probability of traverse equal to zero; however, the algorithm, as
may be easily verified, will place both units at node 2 with probability one giving the evader a successful
traverse probability of 1/16. For one unit, however, the algorithm will optimally place it at node 2
giving the evader a successful traverse probability of 1/4.
N
4 8
FIGURE 2. Sample problem network with node 1 as source and node 9 as sink.
INTERCEPTION IN A NETWORK 215
V. EXAMPLE
For illustrative purposes, the algorithm is used to find the optimal placement of forces for the
network of Fig. 2, with the results of each iteration. Ten intercepting units are available, and all arcs
and nodes have identical probabilities associated with them. Thus, for each node, i, and each arc
(i, /), we have:
p(
p(
p(
p(
p(
p(
P(
P(
P(
) , = 0.6666 p(i,j)i = 0.6666
) 2 = 0.8332 p(i,j) 2 = 0.8332
) 3 = 0.8748 p(i,» 3 = 0.8748
) 4 = 0.8852 p(i,y) 4 = 0.8852
) 5 = 0.8878 p(i,y) 5 = 0.8878
)« = 0.8884 p(i,y) fi = 0.8884
) 7 = 0.8886 p(i, j) 7 =0.8886
) H = 0.8887 p(i, y)« = 0.8887
) 8 = 0.8888 p(i,» » = 0.8888
p(0 ,0 = 0.8889 p(i,j) ,0 = 0.8889
The problem terminates in eight iterations. The results of these iterations follow:
Iteration 1: All arcs and nodes have capacities of 1.5. A maximum flow is 1 unit along path 1,2,5,6,9.
A minimum cut set is node 1 and its value is 1.5. Increasing intercepting units at node 1 to 1 gives a
maximum path value of 0.33333.
Iteration 2: The capacity of node 1 is increased to 2.0; all other capacities are unchanged from
iteration 1. A maximal flow is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 5 and its value
is 1.5. Increasing intercepting units at node 5 to 1 gives a maximum path value of 0.11111.
Iteration 3: The capacity of node 5 is increased to 2.0; all others are unchanged. A maximal flow
is 1.5 units along path 1,2,5,6,9. A minimum cut set is node 9 and its value is 1.5. Increasing intercepting
units at node 9 to 1 gives a maximum path value of 0.03704.
Iteration 4: The capacity of node 9 is increased to 2.0; all others are unchanged. A maximal flow
is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 1 and
its value is 2.0. Increasing intercepting units at node 1 to 2 gives a maximum path value of 0.01853.
Iteration 5: The capacity of node 1 is increased to 4.0; all others are unchanged. A maximum flow
is 1.5 units along path 1,2,5,6,9 and 0.5 unit along path 1,3,5,8,7,9. A minimum cut set is node 5 and
its value is 2.0. Increasing intercepting units at node 5 to 2 gives a maximum path value of 0.00927.
Iteration 6: The capacity of node 5 is increased to 4.0; all others are unchanged. A maximal flow
is 1.5 units along path 1,2,5,6,9 and 0.5 unit on path 1,3,5,8,7,9. A minimum cut set is node 9 and its
value is 2.0. Increasing intercepting units at node 9 to 2 gives a maximum path value of 0.00464.
Iteration 7: The capacity of node 9 increased to 4.0; all others are unchanged. A maximal flow is
1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set is arcs (5,6) and
(5,8) and its value is 3.0. Increasing intercepting units to 1 at arcs (5,6) and (5,8) gives a maximum path
value of 0.00155.
Iteration 8: The capacities of arcs (5,6) and (5,8) are increased to 2.0; all others are unchanged.
A maximal flow is 1.5 units along path 1,2,5,6,9 and 1.5 units along path 1,3,5,8,7,9. A minimum cut set
is nodes 6 and 8 and its value is 3.0. Increasing intercepting units to 1 at nodes 6 and 8 gives a maximum
path value of 0.00052. Algorithm terminates as all 10 units have been assigned.
216 R D  WOLLMER
The final assignment is 2.0 intercepting units each at nodes 1,5, and 9; and 1.0 intercepting unit
each at nodes 6 and 8 and arcs (5,6) and (5,8).
REFERENCES
[1] Durbin, E. P., "An Interdiction Model of Highway Transportation," The RAND Corporation, RM4945PR (May 1966).
[2] Ford, L. R., Jr. and D. R. Fulkerson, "Flows In Networks," The RAND Corporation, R375 (Dec. 1960).
[3] Wollmer, R. D., "Maximizing Flow Through A Network With Node and Arc Capacities," Transportation Science 2, 213232
(1968).
[4] Wollmer, R. D., "Removing Arcs From A Network," Operations Research 12, 934940 (1964).
[5] Wollmer, R. D., "Stochastic Sensitivity Analysis of Maximum Flow and Shortest Route Networks," Management Science
9, 551564 (1968).
[6] Wollmer, R. D., "Algorithms For Targeting Strikes in a Lines of Communication (LOC) Network", Operations Research, 18,
497515 (1970).
THE FIXED CHARGE PROBLEM*
David I. Steinberg
Washington University
St. Louis, Missouri
ABSTRACT
The fixed charge problem is a nonlinear programming problem of practical interest in
business and industry. Yet, until now no computationally feasible exact method of solution
for large problems had been developed.
In this paper an exact algorithm is presented which is computationally feasible for large
problems. The algorithm is based upon a branch and bound approach, with the additional
feature that the amount of computer storage required remains constant throughout (for a
problem of any given size).
Also presented are three suboptimal heuristic algorithms which are of interest because,
although they do not guarantee that the true optimal solution will be found, they usually
yield very good solutions and are extremely rapid techniques.
Computational results are described for several of the heuristic methods and for the branch
and bound algorithm.
1. INTRODUCTION
The fixed charge problem is perhaps one of the most intriguing problems of mathematical pro
gramming. Its structure is almost identical with that of a linear programming problem; yet, the one
difference — the existence of the fixed charges in its objective function — has prevented the development
of any extensive theory concerning its solution.
There are two known exact methods of solution: total enumeration of the extreme points of the
convex set of feasible solutions, and a mixed integercontinuous variable formulation. However, total
enumeration is impractical for large problems, and, at this time, methods for solving mixed integer
continuous variable problems involving a large number of integer variables are not effective enough
to be of much use, either. Therefore, heuristic methods for obtaining a "good" (although not necessarily
optimal) solution have been developed. Some of these methods will be discussed in Section 2.
Also presented is an exact algorithm, based on a branch and bound approach, which has been
developed by the author, along with several new heuristic methods. Computational results are given
in Section 4.
1.1 Formulation of the Fixed Charge Problem
The fixed charge problem is a nonlinear programming problem of the form:
(1) Minimize z = cx + dy
subject to: Ax — b,
(2) x ^ 0,
This research was performed under AEC Research Contract No. AT(111)1493. Submitted originally as Report No.
COO149315.
217
218 D. I. STEINBERG
where A is an m X n matrix of real numbers,
b is an mvector in an mdimensional space, E m ,
c and d are nvectors in E",
x is an rcvector of unknowns, and
y is an nvector such that
fO, if*, =
^ ^={i, if,;> 
It is assumed that all dj ^ 0. The quantities dj are the fixed charges which are incurred if the corre
sponding Xj are (strictly) positive. If all dj = 0, then (1) and (2) reduce to a linear programming problem.
It is known [7] that the optimal solution to the fixed charge problem occurs at an extreme point
of the convex set of feasible solutions (2). Since the total number of extreme points is finite, the problem
could be solved by evaluating the objective function (1) at every extreme point and taking the minimum
of these quantities. However, as is the case with linear programming, the number of extreme points
increases combinatorially with m and n. Hence, it is desirable to develop a method that requires that
only a subset of the extreme points be examined in order to determine the optimal solution. Two such
methods are the mixed integercontinuous variable integer programming formulation and the author's
branch and bound algorithm.
1.2 A Mixed IntegerContinuous Variable Approach
Hadley [7] presents the following formulation of the fixed charge problem and suggests that it be
solved by use of Gomory's mixed integercontinuous variable algorithm.
Minimize z = cx + dy,
subject to: Ax = b
XjUjyj^O, 7=1, . . ., n
0<ft<l, 71 n
yj integers
x^O,
where the Uj represent known upper bounds for the xj.
If A ismXn, then the above formulation has m + 2n constraints. It should be noted, however, that
the upper bound constraints on the yj need not be expressed explicitly, since no optimal solution would
have yj > 1. Thus, the actual basis size required to use Gomory's algorithm is m + n. However, compu
tational experience with Gomory's algorithm has been very disappointing, even for some fairly small
problems [7]; therefore, this approach is not widely used to solve fixed charge problems.
2. SOME HEURISTIC METHODS
The fact that the optimal solution to the fixed charge problem lies at an extreme point of the
convex set (2) has led to the development of some heuristic approaches which use adjacent extreme
point techniques to obtain good solutions. All of the methods discussed below are basically adjacent
extreme point methods.
There are also several heuristic methods for solving the fixed charge transportation problem.
These are presented in [2] and [3], but shall not be discussed here.
FIXED CHARGE PROBLEM . 219
The following notations are offered prior to presenting the heuristic methods:
Xb~ vector of basis variables,
c B = vector of prices for variables in the basis,
aj=yth column of A,
B= basis matrix,
yj — B~ l a,j, and
From the theory of linear programming [6, 7] it is known that once a vector a, is chosen to enter
the basis, the vector a B r which leaves the basis is determined by
XBr
(4) — = Min { XBklykj , for y kj > 0} .
Yrj k
If all dk = (i.e., the problem is all linear) the greatest decrease in the objective function at each
iteration is obtained if the vector to enter the basis, a,, is chosen so that
X t
(cj — Zj) =Min  — {c k — z k ), for y rk > [
k [yrk J
When this expression is nonnegative the optimal solution has been reached. For the fixed charge
problem (one or more dj ^ 0) criterion (5) must be modified to account for the fixed charges; the greatest
decrease in z is obtained if a, is chosen so that
(6) ej = Min {e fc , for y rk > 0},
k
where
(7) e k = — (ck~ z k ) + (d Br — d k )
Vrk
Equation (7) must be modified slightly in the presence of degeneracy, in an obvious way. (It should
be noted that the variable which leaves the basis, x Br , will, in general, be different for each,/). Because
the objective function is concave, it is not true that when all e* 2s a global minimum has been reached.
This stage may merely be a local minimum. That is, no adjacent extreme point yields an improved
value of r, but it is still possible that some other extreme point of the convex set will yield a better value
of 2. The methods discussed below differ in their approach to the problem of trying to obtain a better
solution once a local minimum has been reached. The first three methods are those of the author.
2.1 Heuristic One
1. By using a Simplex Algorithm, vectors are removed from and entered into the basis according
to (4) and (6), respectively, until all e, 2= 0. Let the current solution be denoted by xo, with corresponding
value of the objective function zo. Store xq. Define two parameters:
ao = maximum number of consecutive iterations allowed to find a new extreme point with a value
of 2 smaller than zo.
/3o = maximum number of consecutive iterations allowed in which the objective function does not
decrease.
Also, define two "counters," ai and /3i, which indicate the current number of consecutive itera
tions for which zo has not been improved, and in which the objective function does not decrease, respec
tively. Set ai = /3i = 0.
220 D. I. STEINBERG
2. Find the e, which yields the smallest increase in z and insert a, into the basis. Set ai = /3i = l.
Proceed to Step 3.
3. Find the e, which yields the largest decrease in z and insert a, into the basis. Proceed to Step 4.
4. Compare the new z, z\, with zo:
a) If Zi < Zo, return to Step 1.
b) If z\ = Zo, but Xi 9^ Xo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate.
c) If Z\ = Zo and *i = x , set «i = ai + 1 and /3i = /3i + 1. If ai =£ « and /3i =£ /3 go to Step 5. Otherwise
terminate.
d) If z\ > zo, set ai = ai + 1. If ai =£ ao return to Step 3. Otherwise terminate.
5. Perform /3i consecutive interations in each of which the vector which yields the largest increase
in z is inserted into the basis. Return to Step 3.
This algorithm attempts to find better solutions by moving away from the current local minimum
in as few iterations as possible. The first such attempt (Step 2) actually increases the objective function
as little as possible; if this attempt fails, the objective function is increased as much as possible (Step
5) with the hope that upon once again decreasing (Step 3) a better local minimum will be found. The
efficiency of this method, as well as its accuracy, depends on the parameters ao and /3o. In general, ao
should probably be somewhat larger than /3o. Since computational experience [4] has shown that the
average number of iterations required to solve a linear programming problem with m constraints is
about 3/2m, a reasonable range for a might be a e[m/2, m]. As the computational results of Section 4
indicate, much smaller values may also be used fairly successfully. Of course, economic considerations
may also play a part in the determination of ao and j8o.
2.2 Heuristic Two
1. Same as in previous heuristic, except that the entire Simplex tableau corresponding to xo is
saved. Define a as before, set /3 = n — m, ai = l, /3i=l, and proceed to Step 2 (assuming again that
all ej 3= at xo).
2. Beginning with the tableau corresponding to xo, find the e, which yields the 03i)th smallest
increase in z and insert a, into the basis. Set /3i = /3i + 1, and if /3i =S (J3o + 1) proceed to Step 3. Other
wise terminate.
3. Find the ej which yields the largest decrease in z and insert a, into the basis. Set ai = ai + l
and if ai =£ ao, proceed to Step 4. If ai > a , proceed to Step 2.
4. Compare z\ with zo:
a) If z\ > zo, return to Step 3.
b) If z\ — zo, but jci t^ xo, return to Step 3.
c) If Zi = zo and x\~ xo, set ai = 1 and return to Step 2.
d) If z\ < zo, return to Step 1.
2.3 Heuristic Three
Heuristic Three employs a Monte Carlo approach to find an extreme point with a better value of
z than xo (where again at xo all ej 3= 0). In this method, extreme points within a certain "distance"
from Xo are randomly examined. "Distance" here shall mean the number of basis changes required to
move from one extreme point to another. Let 8 = maximum distance from xo allowed and let ao again
be the maximum number of iterations allowed to find an improvement in z.
The procedure is as follows:
1. Insert a.j into the basis, where ej is the maximum decrease in z. Continue in this manner until
FIXED CHARGE PROBLEM 221
all e, S 5 0. Denote the current extreme point by x . Set a.\ =0 and proceed to Step 2.
2. Randomly select (according to the procedure described below) a vector within a distance
8 of Jto. Denote this vector by X\. Set a t = a t + 1 and proceed to Step 3.
3. Compare z\ with zo:
a) If z\ > zo and «i 3= a , return to Step 2. If z x > zo and «i > a , terminate.
b) If z\ ^ zo, return to Step 1, replacing xo and zo by *i and zi, respectively.
In order to complete the description of the above algorithm, the procedure for the random selection
of extreme points must be described. For simplicity, the vectors selected to enter the basis are chosen
one at a time (and the corresponding vector to leave is determined by (4)). The basis transformation
is then performed before the next vector to enter the basis is again chosen randomly. Since there are
(n — m) nonbasis vectors, the random selection of a vector to enter the basis is performed by generating
a random number from a uniform random distribution on the interval (0, 1), multiplying it by (n — m),
and rounding the result to the next highest integer, say k. Then the A:th nonbasic vector is selected to
enter the basis. The process is repeated 8 times, and the resulting extreme point is the first randomly
generated extreme point within a distance 8 of Xo to be used in Step 2.
An approach very similar to the above is proposed by D. Denzler [5]; the basic difference is that
Denzler's method does not restrict its random search for an extreme point within a certain distance
of the current best solution. Many iterations are thus wasted. Intuitively, it would seem that extreme
points close to one which has a known good value of z should also yield fairly good values of z, and so
it would appear to be a better procedure to restrict the search for better solutions to extreme points
fairly close to Xo rather than attempt to search randomly the entire set of extreme points.
2.4 The CooperDrebes and CooperOlson Heuristics
Cooper and Drebes [3] have developed two other adjacent extreme point heuristics, which differ
from the above three methods mainly in that in the CooperDrebes methods the objective function is
modified at certain stages in the algorithms; in particular, the linear costs cj,j= 1, . . ., n, are modi
fied according to
(8) c j = c j + ^,x j >0
x j
Cj=Cj, xj = 0,
where "*/' implies the current value of variable Xj. At other stages in their heuristics, Cooper and
Drebes remove from the basis the vector with the largest fixed charge (of the current basic variables)
and insert into the basis the vector with the smallest fixed charge (of the current nonbasic variables).
The two methods differ only in their decision rules for employing the various techniques discussed
above.
A modification of the CooperDrebes heuristics has been developed by Cooper and Olson [4]. The
modification consists of randomly generating an extreme point a distance not more than 8 (where
distance is defined as in Section 2.3) from the current local minimum, and then the CooperDrebes
heuristic is again employed in an attempt to improve the solution. The random extreme point generation
is called a "perturbation." The parameter 8 is allowed to vary from about m/3 to m/2. If the search
following the last previous perturbation yields a local minimum worse than the best obtained so far,
8 is reduced; if the objective function value is the same as the best, 8 is increased; if a better value is
obtained, 8 is not changed.
222 D. I. STEINBERG
The computational experience reported on these methods [3, 4] indicates that they will yield the
optimal solution a high percentage of the time and, moreover, can be expected to provide at least a
nearoptimal solution.
2.5 A Separable Programming Approximation Technique
Another technique which yields a local minimum is to formulate the fixed charge problem as a
separable programming problem [7].
Letting Uj represent known upper bound values of Xj, consider the following representation of
the objective function (1):
(9) Z = J£ [djXy +(dj + CjUj)k 2j ] ,
where
d°) xj = \ 2J Uj
and
(11) \0j+klj+\2j=l
and where no more than two \kj can be positive for a given j, and these can only be positive if they are
adjacent. The complete formulation becomes:
N
Minimize z= ]£ [djkij+ (dj + CjUj)\ 2 j],
subject to
^ aijUjkij — bu i = 1 ,
m
AfcjSsO, k = 0, 1, 2 and j=l, . . ., n.
When this problem is solved using a restricted basis entry Simplex method [7], the optimal solution
will be a local minimum for the fixed charge problem. There is no evidence, however, that this procedure
is any better than solving the fixed charge problem by using linear programming and ignoring the fixed
charges. Moreover, no attempts are made to improve the local minimum obtained.
3. A BRANCH AND BOUND ALGORITHM
In attempting to develop an exact solution method for the fixed charge problem, it would be natural
to begin by investigating the possibility of developing an adjacent extreme point method. However,
since the objective function is not always monotonic (that is, there may not always exist a monotonically
decreasing adjacent extreme point path to the optimal solution) some other approach is called for. In
this section a branch and bound approach is presented. (A general knowledge of the branch and bound
philosophy will be assumed in the following discussion. The reader is referred to [1] for a general dis
cussion of branch and bound methods.)
The algorithm finds the global minimum solution without having to examine all basic feasible
solutions of (2). In the next two sections, the general approach will be presented and the computational
details will be given in Section 3.3.
FIXED CHARGE PROBLEM
223
3.1 The Tree of Solutions
Suppose the following tree is generated:
The nodes are numbered from left to right and from top to bottom, so that each horizontal "line"
of nodes may be thought of as a level; that is, the first level consists of only node 1, the second level
consists of nodes 2 and 3, the third level consists of nodes 4, 5, 6, and 7, etc. The constraints (2) will
be required to be satisfied at every node. In addition, each node (except node 1) will be required to
satisfy some additional constraints, of the form "xj=0" or "xj > 0." Since only basic feasible (extreme
point) solutions are of interest, the latter constraints actually require variable Xj to appear in the basic
solution, at a positive level. At each node (except node 1) one additional constraint is imposed, so that
the solution obtained at any given node k must satisfy this new constraint plus the constraints imposed
at every other node in the path from node 1 to node k. Thus, a node at level L must satisfy {L— 1)
additional constraints.
If a path is terminated when the additional constraints imposed by the node are such that:
1. There are m constraints of the form "xj > 0", or
2. There are (n — m) constraints of the form "xj = 0", or
3. No feasible solution exists when the current additional constraints are imposed,
then it is clear that the collection of terminal nodes will contain all the basic feasible solutions to (2).
That is, for terminal case 1, m ay's are required to be positive and hence a basis is uniquely determined.
For case 2, (n — m) x/s must be zero, and hence the remaining xfs uniquely determine a basis; however,
this basis may not yield a feasible solution, in which case terminal case 3 applies.
Now, of course, the value of the objective function z can be computed for each extreme point and
the minimal value obtained; however, it is desirable to reduce as much as possible the number of
extreme points which must be examined.
3.2 The Algorithm
The following notations will be introduced:
S* — {j\ x j must be positive at node k},
Sj = {j\xj must be zero at node k},
' z k — a lower bound of 2 at node k, and
z* = the global optimal value of z.
224 D  ! STEINBERG
The constraint set at node k thus becomes:
/ Ax = b
(13) T k =
x ^
Xj > 0, >eSf
Xj = 0, jeS k
x is an extreme point.
Suppose that an upper bound estimate of z*, z, is obtained. Then a path can also be terminated
at node k if
(14) z k =s z,
since the basic feasible solution corresponding to z has as good a solution as any of the solutions which
could be obtained by continuing on the current path beyond node k.
This procedure will cut down on the number of basic feasible solutions that need to be examined;
its efficiency depends heavily on how good the bounds z k are and on how close z is to the true optimal
value of z.
Any heuristic method can be used to obtain an initial z\ thereafter, every time a path is terminated
because a unique basic feasible solution has been obtained, the corresponding value of z, z, can be
calculated and compared with the current z; if z is smaller than z, then z and its corresponding x
replace z and its corresponding x as the current upper bound estimate of the optimal solution. A method
for calculating z k is given in the next section.
The branch and bound algorithm can be summarized as follows:
1. An initial upper bound, z, of the optimal value of the objective function is obtained by means
of a heuristic. This bound and its corresponding solution vector x are stored.
2. The paths are then generated so that in branching from one node to another, the node at which
the additional constraint is of the form "xj > 0" is selected as the next node in the path, until one of
five path terminating rules occurs. The first three such rules are given on page 223. Rules 1 and 2 in
dicate that a unique solution has been reached (although it may not be feasible). Rule 3 pertains to
infeasibility; that is, no basic feasible solution exists for the constraints (13). Rule 4 is the initial bound
ing procedure: A path is terminated if the lower bound z k of the objective function is as large as the
current upper bound of z*; that is, if (14) holds, then the basic feasible solution corresponding to z
has as good a solution as any of the solutions which could be obtained by continuing on the current
path. Hence, the path is terminated. Rule 5 pertains to the fixed charges themselves. When the only
fixed charges whose subscripts are not in either Si or S k are all zero, then as will be seen in the next
section, the lower bound z k is the actual minimum value of the objective function (1) subject to the
constraints (13). Hence, the path can be terminated (and if z k < z, then z k and its corresponding solution
vector become the new upper bound).
3. When the last path has been terminated, the current lower bound z and its corresponding x will
be the optimal solution.
3.3 Computation of the Lower Bounds
In order to obtain the lower bound z k some initial computations are made. In particular, it will
be useful to know the level of degeneracy that exists in (2); that is, how many basic variables can
appear in a basic solution at a zero level. To obtain this information, the following n linear programming
problems are solved:
FIXED CHARGE PROBLEM . 225
Minimize Zj = Xj,
subject to: x is an extreme point of the convex set (2) and Xj is in the basis.
The method for solving these restricted linear programming problems will be presented later in this
section.
If the optimal solutions to (15) are denoted by xf, respectively, then if all xf > 0, no degeneracy
exists; if any x*—0, then the number of such variables determines the maximum level of degeneracy
which could occur.
It shall be assumed that:
(16) di^d2^ . . . ^d„.
In practice (16) is really no restriction, since the vectors may only have to be reordered so that (16)
holds.
Let S$ = {j\j4S* US*},
P—{j\xj can be in the basis at a zero level},
c?=pn (Sfusf),
N k = the number of elements in S k : , and
N k i = the number of elements in Q.
Thus S k is the set of subscripts whose corresponding variables do not appear as an additional
constraint at node k; the set Q consists of those subscripts whose corresponding variables could appear
in the basis at a zero level at node k.
Consider the following problem:
(17) Minimize z\~cx
subject to: xeT k .
The desired optimal solution to problem (17), for any node k, lies at an extreme point of the constraints
(2) and also satisfies the additional constraints imposed at node k. Denote this optimal solution by
z k and let
(18) 2§ = 2 dj.
Further, let z denote the value of the sum of the last (m — N k — N k i ) fixed charges. Define
(19) z k = z k + z k + z k .
In order for z k to be a legitimate lower bound for z*, it is sufficient to show that for any solution xeT k
(i.e., any solution x which satisfies the constraints at node k),
(20) cx + dy^z k .
But,
cx + dy=cx+ V g?j + Y djyj,
(21) ;eSf Ui
since for./€Sf, yj= 1 and for^'eS^', yj— 0. By definition, ex 3 s z k .
(22)
226 D. I. STEINBERG
Since the maximum level of degeneracy possible at node k is N£, there must be at least (m — N£) fixed
charges whose yj = 1 in any basic solution at node k; but,N k of these fixed charges are determined by S^
(and are included in z k ). Thus, there are at least (m — N k — Nq) remaining fixed charges whose yj must
be unity. Since z k is defined to be the sum of the (tii — N* — Nq) smallest remaining fixed charges, it
must be true that
(23) 2 d M * 2 "
its*
Addition of the three relations (18), (22), and (23) yields the desired result (20). Hence z k as defined by
Eq. (19) is a legitimate lower bound for the optimal solution which could be attained at node k.
In order to complete the description of the algorithm it is only necessary to describe the procedure
for solving the linear programming problems (15) and (17). Since (15) is merely a special case of (17)
only the latter is considered below.
Consider the following problem:
(24) Minimize Z\ = ex
(25) subject to: xeT k .
The set of extreme points which satisfy the constraints (25) is a connected set; that is, an adjacent ex
treme point path exists between any two points in the set. Suppose an initial basic feasible solution to (25)
is known, and that a Simplex algorithm is used to improve the objective function. Then, in a finite number
of steps, a stage will be reached (call it Stage One) in which the only vectors which could be entered into
the basis to yield a decrease in the objective function would require one of the vectors whose corresponding
jeS k to leave the basis. (It is assumed throughout the discussion that none of the vectors whose corre
sponding jeS k are considered for entry into the basis.) Let the current solution at this stage be denoted
by x, with corresponding i\ = ex. It is possible that some other extreme point £o exists such that zo < i\ ,
but that x and xo are not adjacent extreme points. In order to determine whether or not 5c is indeed an
optimal solution to problem (17), the following constraint is appended to (25):
(26) cx + x n + l = z l .
The basis size is now increased to m + l, and the (m+ l)st basic variable for the current solution is
x n + i = 0. Now, however, it is possible that one or more variables may be brought into the basis without
removing any of the variables whose corresponding jeS k . If $ is not an optimal solution to (17), then a
solution exists for which x n +i > 0. In order to determine whether such a solution does in fact exist,
the following test for optimality is employed:
Let R k denote the set of variables which are not in the basis for x (and whose j4S k ) . Then, for each
xi in R k , first insert xi into the basis at a zero level, replacing x n + u and then enumerate the extreme
points (of basis size m + 1 ) of the convex set determined by (25) and (26), by solving
Maximize z = x,,
subject to: Ax = b
CX + Xn+l = Z\
(27) Xi — x„ +2 = xt
x ^
FIXED CHARGE PROBLEM • 227
xj > 0, jeS$
Xj = 0, jeS k
x is an extreme point,
where x, is the current value of x,, initially zero. At each step in solving (27) the minimum increase in
Xi is taken. In this manner every extreme point of (25) and (27) which contains x, will be examined. If
at any stage x n +i > 0, the procedure can be terminated, and by beginning with this new extreme point,
the objective function (24) can be reduced in the usual way until Stage One is again reached. If, on the
other hand, the optimal solution to (27) for each x\ is zero, or no extreme points for which x n + \ > have
been found after solving (27) for all XitR k , then x is the optimal solution to problem (17).
This procedure will not, in general, be as lengthy as it might at first appear. If x is the unique
optimal solution to problem (17), then the problems (27) will all have optimal values of zero, and one
iteration per problem will indicate this. Moreover, if x is not in fact optimal, is should not require very
many iterations to find a solution for which x n +i > 0.
It should also be noted that the true minimal solution to problem (17) is required only when
i\ > (z — z k — z k ), where z is the upper bound estimate, since if z.\ =£ z the path will not be terminated
and decreasing z\ further is of no interest. Hence the above test for optimality for problem (17) need
be applied at any node only when i\ > z.
The question still remains: How is an initial basic feasible solution to (17) obtained? The procedure
is as follows:
Let W\ = {xj\jeS k , x } not in current basis}
W<i = {xj\J€S%, Xj > in current basis} .
Thus, those XjeWi satisfy the constraints imposed at node k, and those XjeWt do not; in order to remove
those XjtWt from the basis, the following linear programming problems may be solved, for each xfiW 2 :
Minimize z=Xj,
(28) subject to: Ax= b
x^O.
When z = 0, either xj has been removed from the basis, or it is in the basis at a zero level and can be
removed; if in an optimal solution z > 0, no feasible solution exists with Xj not in the basis.
It should be noted that once a variable XjeW 2 has been removed from the basis, then XjeW\\ from
the theory of linear programming it is known that if it is desired that some variables not appear in a
basic solution, they may merely be excluded from consideration for entry into the basis; thus, pro
ceeding to remove the variables XjeW 2 from the basis while never considering the variables XjeW\ for
entry into the basis will yield in a finite number of steps a basic solution in which W\ = S^ and W 2 is
empty.
Now, the variables XjtS'l must be inserted into the basis:
Let Fi = {xj\jeS\ and xj > 0}
Y 2 = { Xj \J€S k 1 and Xj = 0}.
228
D. I. STEINBERG
For those variables XjeY 2 the following linear programming problems are solved:
Maximize z — Xj
subject to: Ax = b
(29) x ^
XjeYi remain in the basis
at a positive level
x is an extreme point.
Since the XjeY x are already in the basis at a positive level, the problems (29) can be solved by the
method discussed earlier in this section. If the optimal solution to any of the above problems is zero,
then no feasible solution to problem (29) exists with that variable in the basis at a positive level; if any
solution exists for which z > 0, then variable xj has successfully entered the basis (it is not necessary
actually to optimize the objective function). At this point, Xj will now be in Yi and the next variable
XjeY 2 can be inserted in the basis via problem (29). In a finite number of steps, a stage will be reached
at which Y\ = S\ and Y 2 is empty. If none of the variables XjeS% were considered for entry into the basis
during the procedure of solving the problems (29), then an initial basic feasible solution to problem
(17) has been obtained.
It should be noted that at no stage of the branch and bound procedure will more than one variable
not satisfy the additional constraints imposed by any node k. Thus, the amount of actual computation
necessary to obtain the initial solution is much less than one might at first suspect.
Moreover, the amount of computer storage required by the algorithm does not build up as the
number of nodes examined increases, but remains constant throughout the procedure (for a problem
of any given size). This advantage does not exist in most other branch and bound algorithms.
4. COMPUTATIONAL RESULTS
The primary difficulty in finding fairly large problems to test heuristic methods is that the optimal
solution must be known. Cooper and Drebes [3] have randomly generated 290 (5 X 10) problems with
the following properties:
(1) ay*£20,
(2) l=£c jS £20,
(3) 1 *£ dj *£ 999,
(4) The density of A is 50 percent.
The optimal solutions to these problems were obtained by total enumeration.
In order to take advantage of these problems, larger problems were constructed by forming the
A matrix as follows:
""/I, • • •
A 2 ■ • •
A =
A k
FIXED CHARGE PROBLEM 229
where the Ai are (5 X 10) matrices. The resulting problem is of size (5k X 10k) and its corresponding
solution is known. Hence problems of varying size can be formulated to test the heuristic methods.
Because of roundoff difficulties, some of the test problems are not included in the results. These are
noted under the proper table.
Since the heuristics tested are all adjacent extreme point methods, each iteration corresponds to a
regular linear programming transformation. Thus, the primary means of comparison among different
heuristic approaches is the number of iterations before the best obtained solution is found. Since these
methods are not exact, the "best obtained solution" will not always be the optimal solution, and in
some cases will be a function of the parameters related to the particular method.
The approximate computer times (for an IBM 360/50) per iteration for the four groups of problems
for both Heuristic One and Heuristic Two were (in seconds): 0.10, 0.18, 0.60, and 2.4.
The Cooper and Drebes paper [3] reports computational experience with two heuristics: both
heuristics were tested on a set of 253(5 X 10) problems; 240 were solved exactly by the first method
in an average of 75 total iterations per problem, and 245 problems were solved exactly by the second
method in an average of 100 total iterations per problem. The three (5 X 10) problems which neither
heuristic solved exactly were all less than 10 percent from optimality. No information on the number
of iterations required to achieve the best solution is given.
On a set of 70(15 X 30) problems, the two methods were combined; the resulting heuristic solved
63 problems exactly with an average of 1200 iterations per problem. The nonoptimal problems were
less than 3 percent from optimality.
4.1 Heuristic One
Problems of size (5 X 10) ,( 15 X 30) , (30 X 60) , and (50 X 100) were used in the testing of Heu
ristic One. Of the 250(5 X 10) problems, 235 were solved exactly with olq = 2 and /3o = 2. The average
number of iterations required to achieve the final solution was 3 and the average total number of
iterations per problem was 6.8. Table 1 contains a summary of the 15 nonoptimal problems. When ao
was increased to 3 and /3o to 4, problems 2, 94, and 128 were solved exactly. No further increase in
ao and /3o yielded any improvements in the 12 remaining suboptimal problems.
Ninety (15 X 30) problems were attempted, and with a = 2, /3 — 2, the optimal solutions to 75
problems were obtained, with an average of about 5 iterations before the best solution was obtained
and an average of 12 total iterations per problem. No further improvement resulted when ao was
increased to 10; the average total number of iterations increased to 16. /3o was increased again to 15,
with an improvement in only one of the suboptimal solutionsproblem 43. The optimal solution is
Z* = 8317.7 and the best previously obtained solution was Z = 9172.4; with /3 = 15, the solutions
8964.7, 8593.0, and 8525.8 were obtained. The average total number of iterations was 23. Table 2
contains a summary of the suboptimal (15 X 30) problems.
230
D. I. STEINBERG
Table 1. Heuristic One(5 X 70) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number"
solution
solution
from optimal
2
2830.4
"3317.6
14.7
57
4015.8
4017.9
0.05
70
3363.0
3468.3
3.1
86
2359.4
2380.7
0.9
93
2351.7
2436.6
3.6
94
2373.0
"2899.1
1.1
99
2418.0
2473.0
2.3
128
3459.0
» 3470.6
0.3
131
1543.0
1766.1
14.5
137
3593.1
3885.1
8.1
181
2826.2
3183.0
12.6
188
2417.1
2440.1
0.95
235
2080.9
2155.9
3.6
267
3072.2
3170.9
3.2
279
2323.4
2533.7
9.1
"Problems not attempted: 29, 36, 51, 55, 66, 67, 73, 75, 82, 88, 98, 117, 119, 120, 126, 130, 134, 141, 147, 158, 162, 164, 166,
176, 185, 201, 204, 206, 211, 215, 219, 226, 236, 237, 238, 248, 250, 260, 264, 266, 275.
b Exact solutions obtained with a = 3 and /3o = 4.
Table 2. Heuristic One(15X30) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number"
solution
solution
from optimal
20
7676.4
7726.9
0.67
29
7522.6
7544.0
0.28
30
9452.9
9458.2
0.06
31
9201.7
9286.5
0.89
32
6666.4
7192.6
7.89
33
8738.3
8793.3
0.63
40
11,117.3
11,365.9
2.23
43
8317.7
8525.8
2.50
45
6785.6
6938.3
2.25
46
8517.7
8809.7
3.43
53
6497.6
6647.7
2.31
61
7651.0
8007.7
4.66
63
7836.0
7859.6
0.29
85
6307.2
6635.2
5.20
93
8592.7
9047.5
5.29
" Problems not attempted: 1, 44, 64, 68, 79, 89.
FIXED CHARGE PROBLEM
231
Thirtysix (30x60) problems were attempted and 26 were solved optimally with a — 6 and /3 = 25.
The average number of iterations before the best solution was obtained was about 11 iterations and the
total number of iterations per problem averaged about 42. A summary of the suboptimal (30X60) prob
lems is given in Table 3.
TABLE 3. Heuristic One(30 X 60) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number "
solution
solution
from optimal
15
16,975
17,002
0.04
16
15,868
16,479
3.98
17
18,441
18,4%
0.30
23
15,303
15,748
2.91
25
16,886
16,987
0.50
27
13,960
14,110
1.07
31
14,944
15,301
2.39
32
17,350
17,440
0.52
37
14,716
15,181
3.16
43
13,946
14,274
2.35
° Problems not attempted: 1, 5, 10, 20, 22, 24, 34, 40, 45.
Nine(50X100) problems were also tested, with a ~ 10 and /3o = 40; the results are given in Table 4.
Table 4. Heuristic One(50 X 100) Problems
Problem
number
Optimal
solution
Best obtained
solution
Percent
from optimal
No. of iterations
to best obtained
solution
Total No. of
iterations
1
2
5
12
14
15
18
19
20
26,761
28,842
25,545
29,350
24,444
29,380
26,844
25,903
27,822
26,761
28,842
25,545
30,182
25,085
29,481
26,844
26,350
27,899
47
22
21
14
28
13
18
24
18
87
62
61
49
46
48
51
49
48
2.83
2.57
0.34
1.73
0.28
Average No. of Iterations
22.8
55.7
4.2 Heuristic Two
The same sets of problems used to test Heuristic One were employed in the testing of Heuristic
Two. In some cases, again, roundoff difficulties caused some problems to be discarded, and these
problems are again noted in the tables.
232
D. I. STEINBERG
Of the 268 (5X10) problems attempted by Heuristic Two, 255 were solved optimally. A summary
of the suboptimal problems appears in Table 5. The average number of iterations required to achieve
the final solution was 6, and the average total number of iterations was 15.2.
TABLE 5. Heuristic Two(5 X 10) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number "
solution
solition
from optimal
28
2993.4
3025.4
1.1
50
2383.9
2813.4
18.0
71
3058.3
3470.5
13.5
93
2351.7
2436.6
3.6
99
2418.0
2473.0
2.3
108
2414.8
2423.8
0.4
119
3467.6
3560.6
2.7
120
3579.1
3734.5
4.3
131
1543.0
1766.2
14.5
137
3593.1
3885.1
8.1
154
2589.5
2622.7
1.3
245
3574.7
3635.8
1.7
279
2323.4
2533.7
9.1
"Problems not attempted: 23, 27, 30, 52, 61, 67, 76, 81, 94, 98, 110, 128, 149, 158, 165, 202, 210, 224, 225, 239, 249, 254.
Of the 287 (5X10) problems attempted by at least one of Heuristic One and Heuristic Two, 280
were solved optimally; the following 7 problems were not solved optimally by either method: 93, 99,
119, 120, 131, 137, and 279.
Eightyfour (15X30) problems were attempted by Heuristic Two, and 74 were solved optimally.
The average number of iterations required to obtain the final solution was 12, and the average total
number of iterations per problem was 43.
Table 6. Heuristic Two(15 X 30) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number °
solution
solution
from optimal
1
8096.2
8142.6
0.57
31
9201.7
9286.5
0.89
36
6162.5
6171.6
0.15
40
11,117.5
11,365.9
2.23
43
8525.8
8537.4
2.50
44
6807.3
6905.5
1.44
46
8517.7
8809.7
3.43
52
7880.0
8214.1
2.97
79
7765.0
7840.0
2.45
93
8592.7
8803.0
5.29
"Problems not attempted: 16, 18, 21, 22, 26, 27, 33, 37, 61, 72, 80, 83.
Of the 96 (15 X 30) problems attempted by at least one of Heuristic One and Heuristic Two, only
the following 5 problems were not solved optimally: 31, 40, 43, 46, and 93.
,u i^riAKGE PROBLEM
233
TABLE 7. Heuristic Two(30 X 60) Suboptimal Problems
Problem
Optimal
Best obtained
Percent
number
solution
solution
from optimal
1
15,490
15,536
0.30
20
18,407
18,655
1.34
22
15,333
15,535
1.32
23
15,303
15,595
1.91
26
16,848
17,346
2.91
29
15,773
15,889
0.73
34
16,912
17,076
0.97
41
15,601
15,662
0.31
45
14,243
14,267
0.17
47
15,852
16,062
1.32
i
"Problems not attempted: 5, 9, 11, 13, 14, 16, 17, 19, 31, 36, 40, 42, 43.
Thirtysix (30 X 60) problems were attempted by Heuristic Two, and 26 were solved optimally.
The average number of iterations required to obtain the final solution was 25, and the average total
number of iterations per problem was 81.
Of the 46 (30 X 60) problems attempted by at least one of Heuristic One and Heuristic Two, 42
were solved optimally; the following problems were not solved optimally by either method: 1, 23, 34,
and 45.
Heuristic Two was not tested on any larger problems.
4.3 The Branch and Bound Method
Fifteen (5 X 10) problems and 10 (15 X 30) problems were solved by the branch and bound pro
cedure. The results are given in Tables 8 and 9. The total possible number of nodes that might have to be
Table 8. The Branch and Bound Algorithm(5X 10) Problems
Problem
Number
Time
number
of nodes
(minutes)
1
66
0.23
2
24
0.16
3
20
0.16
4
44
0.26
5
20
0.13
6
14
0.11
7
20
0.13
8
74
0.31
9
10
0.09
10
38
0.23
11
34
0.17
12
12
0.10
13
40
0.16
• 14
24
0.13
15
42
0.30
Average
32
0.16
234
D. I. STEINBERG
TABLE 9. The Branch and Bound Algorithm(l 5X30) Problems
Problem
Number
Time
number
of nodes
(minutes)
6
1242
17.9
8
672
9.1
9
3144
47.3
14
1134
18.3
18
308
8.5
29
2294
43.0
34
502
26.4
35
334
4.8
37
904
11.6
38
1546
24.4
Average
1208
21.1
examined is 2" +1 — 1. The average number of nodes that were examined was 32 and 1208, for the two
groups of problems, respectively. Also included in the tables is the computer time required to move
through the nodes; these times are for an IBM 7072. It is felt that these times could be greatly improved
if the computer code were made more efficient. There are many improvements that could be made in
the code, such as rewriting it in machine language and utilizing a more efficient Simplex method (such
as revised Simplex).
The average number of nodes examined for the (5 X 10) problems was 2 5 out of a possible total
of approximately 2 10 nodes; for the (15 X 30) problems the average number of nodes actually examined
was approximately 2 10 out of a possible total of approximately 2 30 . Hence, the ratio of the exponents
of the number of nodes examined to the number of total possible nodes decreased when the number
of constraints was tripled (from 5 to 15). If it is assumed that this ratio continues to decrease linearly,
then the number of nodes that would have to be examined for a problem of size (45 X 90) is about
2 20 nodes.
Hence, it appears that the number of nodes examined grows rapidly with the number of variables,
and so the efficiency of the computer code is extremely important.
5. SUMMARY AND CONCLUSIONS
The computational results of Section 4 clearly indicate that the heuristic algorithms are orders
of magnitude faster than the branch and bound algorithm. Moreover, the heuristics almost always
provide very good, if not optimal, solutions. When two or more heuristics are combined, the likelihood
of finding the optimal solution increases considerably. In fact, when the two heuristics of Cooper and
Drebes [3] are combined with Heuristic One and Heuristic Two, all 290 of the (5 X 10) problems were
solved optimally, only one problem of the 96 (15 X 30) problems (problem 53) was not optimally solved.
The branch and bound algorithm presented herein is a new, exact algorithm for the solution of
the fixed charge problem and appears to be much better than the mixed integercontinuous variable
approach discussed in Section 1.2. Such an exact method of solution is essential for those situations
in which the true optimal solution must be determined.
It is felt that many improvements could be made in the branch and bound algorithm, in addition
to refinements in the computer code used to test the algorithm. It is hoped that a better method of
FIXED CHARGE PROR.EM 235
determining the true maximum level of degeneracy of the convex set of solutions can be found, so that
a better lower bound for the objective function at any given node may be computed. Other refinements
in the handling of the degeneracy problem could probably also be found.
It is suggested that research be conducted on the problem of generalizing the branch and bound
approach described herein, so that it would be applicable to more general separable concave functions.
Another generalization which should be considered is the case in which the fixed charges are dependent
upon more than one variable. For example, a fixed charge dtj could be incurred only if the corresponding
variables xt and Xj are both positive.
6. ACKNOWLEDGMENTS
The author wishes to acknowledge the support of the Washington University Computing Facilities
through NSF Grant G22296 and also the Atomic Energy Commission for partial support of the research
under Research Contract No. AT(111)1493.
REFERENCES
[1] Agin, Norman, "Optimum Seeking with Branch and Bound," Management Science, 13 B 176 185 (Dec. 1966).
[2] Balinski, M. L., "FixedCost Transportation Problems," Nav. Res. Log. Quart. 8 4154 (Jan. 1961).
[3] Cooper, L. and C. Drebes, "An Approximate Solution Method for the Fixed Charge Problem," Nav. Res. Log. Quart. 14,
101113 (Mar. 1967).
[4] Cooper, L. and A. M. Olson, "Random Perturbations and the MIMII Heuristics for the Fixed Charge Problem," Department
of Applied Mathematics and Computer Science, Washington University, Report No. COO14937.
[5] Denzler, D. R., "An Approximative Algorithm for the Fixed Charge Problem," Nav. Res. Log. Quart. 16, 411416 (Sept 1969).
[6] Hadley, G., Linear Programming (Addison Wesley Publishing Company, Inc., Reading, Massachusetts, 1962).
[7] Hadley, G., Nonlinear and Dynamic Programming (Addison Wesley Publishing Company, Inc., Reading, Massachusetts,
1964).
[8] Robers, P. and L. Cooper, "A Study of the Fixed Charge Transportation Problem," Department of Applied Mathematics and
Computer Science, Washington University, Report No. COO 14939.
ON COMPUTING THE EXPECTED DISCOUNTED RETURN IN A
MARKOV CHAIN*
David F. Hitchcock
and
James B. MacQueen
University of California, Los Angeles
ABSTRACT
The discounted return associated with a finite state Markov chain X t , X 2 ■ ■ ■ is given by
g(Xi) + ag(X 2 ) + a 2 g(X 3 ) + . . ., where g(x) represents the immediate return from state
x. Knowing the transition matrix of the chain, it is desired to compute the expected discounted
return (present worth) given the initial state. This type of problem arises in inventory theory,
dynamic programming, and elsewhere.
Usually the solution is approximated by solving the system of linear equations character
izing the expected return. These equations can be solved by a variety of wellknown methods.
This paper describes yet another method, which is a slight modification of the classical itera
tive scheme. The method gives sequences of upper and lower bounds which converge mono
tonely to the solution. Hence, the method is relatively free of error control problems.
Computational experiments were conducted which suggest that for problems with a large
number of states, the method is quite efficient. The amount of computation required to obtain
the solution increases much slower with an increase in the number of states, N, than with
the conventional methods. In fact, computational time is more nearly proportional to /V 2 ,
than to N :i .
I. Introduction
Consider a Markov chain X x , X 2 , . . ., taking values in the finite set S= [1,2, . . ., /V] that has
transition probabilities p(y;x) = Pr[X t+ i — y\Xt = x]. When in state x an amount g{x) is received
immediately. The discounted return associated with passing through the sequence of states
X\, X2, ■ ■ , is
g{X l ) + ag(X i ) + c?g(X 3 ) + . . .,
where a is the discount factor (0^a<l).
The expected discounted return given Xi = x, represented by u(x) (also called the present worth
of state x), is well known to be the unique solution of the equation
(!) u(x) = g(x) + a^ u(y)p(y;x).
An efficient machine algorithm for solving (1) is desired.
The problem of computing the expected discounted return frequently arises in inventory problems,
replacement problems, and elsewhere (c.f. [3], [4]).
If N (the number of states in S) is reasonably small, then (1) can be solved by familiar techniques,
such as Gausian elimination or the simplex method, and is a straight forward problem in machine com
putation. However, if there is a large number of states, the amount of computation, which is generally
*This work was supported by the Western Management Science Institute under a grant by the Ford Foundation, and by
the Office of Naval Research under Contract No. 233(75), Task No. 047003. Computational facilities were provided by the
Western Data Processing Center, University of California, Los Angeles.
237
238
D. F. HITCHCOCK AND J. B. MACQUEEN
proportional to N 3 (cf.[l, p. 215]) using such methods, may be prohibitively expensive. Also, with these
methods error control problems can arise even in problems of moderate size.
This note describes some computational experiments using a slightly modified version of the
classical iterative method for computing u. The method determines sequences of functions {u' n } and
{u'n) on S which converge monotonely to u from below and above, respectively. Then u n = (u' n + u'^)/2
is taken as the estimate of u. The experiments indicate that for "large" problems (say /V ^ 100) , the
amount of computation required to obtain a solution whose percentage error at each state is less than
a preassigned figure, increases much more slowly than /V 3 . For example, in a number of the experiments
with TV = 500, approximately 16(500 2 ) elementary (multiply and add) computations were required to
achieve 0.005percent accuracy. As with the estimate produced by the classical iteration scheme
(see Eq. (5) below) convergence is found to be exponential, but at a rate significantly faster than the
classical method. The monotonely convergent upper and lower bounds, which are easily calculated,
provide a very satisfactory means of controlling error due to roundoff, which may accumulate with
other methods.
The method starts with any function v\ which is at a selected state seS, with the functions t>2,
t>3, . . ., being defined by
(2) v„ + 1 (x) = g(x) + a 2 Vn(y)p(y; x)  (g(s) + a J v n (y)p(y; *) )•
Thus v n (s) = 0, n = 1, 2, . .
(3)
. ... Let
L' n = min x
L" n = max x
g(x) + a ^ v n (y)p(y; x)  v„(x)
y
g(x) + a ^ v n {y)p(y; x)  v„(x)
The sequences {u' n } and {u" n } are then given by
(4) u' n = v„+ (1 a)~ l L' n
u: = v n + (la)'L;
The sequences {u' n } and {u'^\ have the properties:
(i) u' n ss u =£ ul, (ii) u' n s£ u' n + 1 , u" n ^ u'„ +u and (hi) u' n * u, u" n * u.
These properties are obtained directly as special cases of Theorem 2 in [2].
The modification represented by the calculation of u„ is to be compared with the classical iterative
method. The latter is defined by an initial function U\ and the rule
(5)
u„ + i(x) = g(x) + a^ u n (y)p(y; x).
If ui{x) =vi(x) + «i(s), then v n (x) = u n {x) — u n (s), and
«„U) = «„(x) + V2(ia)'(y;+y;'),
y^ = min x g(x) + a^u n (y)p(y; x) — u n {x)
L v
Jl=m<ix x g{x) +a^u n (y)p(y; x) —u„(x)
(6)
where
(7)
COMPUTING RETURN IN MARKOV CHAIN 239
II. Experimental Procedure
To explore the computational efficiency of the suggested method, it was decided to investigate
four main types of transition matrices, with TV equal to 10, 50, 100, 500, or 1,000. The discount factor
was varied over the four values 0.90, 0.95, 0.99, and 0.995 representing the likely range of applications.
In all the experiments, the initial function v\ was identically zero.
The four main types of transition matrices were as follows:
(I) Random. In this type of matrix, the row entries were randomly generated by selecting for each
entry a uniformly distributed random number between and 1 and then normalizing by dividing by
the row total.
(II) Large probability in the left hand column. Here the first column has a specified value po (i.e.,
p(l; x) —po for x= 1,2, . . ., N). The value of po was varied over the three values 0.1, 0.5, and 0.9,
with the remaining entries generated randomly as above (summing to 1 — po, of course).
(III) Large probability on the diagonal. Here all diagonal elements were set equal to po, and the other
entries were random as before. The value of po was varied as above.
(IV) Large probabilities scattered. In this type of matrix, the first row was generated by taking four
probabilities pi=p /2, p 2 = po/4, p 3 = po/8, and p 4 = p /8 and spacing them evenly in the first row.
For example, for /V=100 we have p(l; l)=pi, p(25;l)=p 2 , p(50; l)=p 3 , and p(75;l)=p 4 . The
remaining rows of p(y; x) were obtained from the first row by cyclic permutations, the last entry in
the first row becoming the first entry in the second row, etc.
It was decided arbitrarily to have g(x) vary in multiples of 5 as follows: for n—10, g(l) = 5,
g{2) = 10, . . ., g(10) = 50. For the remaining values of N, all of which were multiples 10, the first
N/10 values of g were 5, the second N/10 values were 10, etc., with the last A710 values being 50. Evi
dence that this method of choosing g, while arbitrary, was not malicious, was obtained by repeating
all the experiments for N= 100 for matrices of Type III (large probabilities on the diagonal) with a
different scheme for choosing g. Each value of g(x) , x— 1, 2, . . ., N, was determined as a uniformly
distributed random number between and 100. Matrices of Type HI were chosen for this additional
experimentation since they were the most difficult of the four in terms of convergence times. The
number of iterations (see below, Table 1) using the randomly determined g was essentially identical
to those obtained with the fixed g, differing by 0, 1 or at most 2 iterations.
III. Results
The main results are presented in Table 1. The entries in the table represent the number of
iterations (the value of n in Eq. (4)) required for the maximum possible percentage error (MPPE) of
the estimate u n to be less than 0.005 or 0.5 percent. Thus, for the indicated value of n, max x \u(x)
— iin(x)  /u(x) is less than 0.00005 or 0.005. The upper and lower entries in each cell in the table
correspond respectively to the 0.005 and 0.5 percent figures. The values of n were actually computed
by allowing the computer program to run until max x {u'n{x) — u' n (x))lu' n {x) =£0.0001, which implies
max x «(x) — u(x)\lu(x) =£ 0.00005. The 0.5percent values were then determined by inspection of
the output which included the MPPE for each iteration.
It should be noted that the MPPE may be considerably larger than the actual percentage error.
As a typical example, consider the matrix of Type III with P = 0.90, /V= 100 and a = 0.90 where the
MPPE after 16 iterations was 6.4 percent. However, ui 6 (x) differed by at most 1.4 percent from u 52 (jc),
and the latter (which corresponds to the last iteration performed) has at least 0.005 percent accuracy.
240 D  F HITCHCOCK AND J. B. MACQUEEN
The MPPE would not be a meaningful criterion of accuracy if u(x) were very near zero for some x.
This possibility was avoided by the choice of strictly positive values of g. The MPPE can also be
rendered meaningless by adding a large positive constant to each value of g, since if u is the solution
with immediate return function g, u+A/(l — a) is the solution with immediate return g + A. Thus by
making A sufficiently large, the MPPE in the approximate and easily computed estimate, 0. — A/(l — a)
can be made arbitrarily small. Again, the choice of g avoids this trivial extreme.
Several general findings were noted. First, the MPPE exhibited a strong tendency to decrease
exponentially with the number of iterations, i.e., as kr n ,0 > r > 1. Moreover k is typically much less
than unity. This was clearly evident on plotting the error against n on semilogarithmic paper. The
maximum percentage error for the classical method defined by equation (5), also decreases exponen
tially, approximately as Ka" where K is near unity. However, r was appreciably less than a, the ratio
r/a being .9 or less in all the experiments. Thus the ratio of percentage errors for the two methods after
30 iterations, say, would be around .05 and perhaps smaller.
A second finding is that the method automatically takes advantage of certain simplifying features
of the transition matrices. Thus in matrices of Type I and II, the distribution of the next value of g
is very similar from state to state, and the sequence g(Xi), g(X 2 ), . . ., is a sequence of similarly
distributed, nearly independent random variables. For that reason, a good estimate for u(x) is simply
the immediate return g( x) plus the mean value of g times a/ (I — a), and a good computational scheme
ought to take advantage of this simplicity. The experimental results show that the method does. Some
what surprisingly, when the rows of the transition matrix are identical, then
U2(x) =u'i(x) =g(x) +a(la) 1 ^g(y)p(y; x)=u(x),
y
and the exact result is obtained immediately (given V\ — 0).
This suggests a general rule, that the less the effect of the initial value of x on the sequence of
ensuing returns, the less the amount of computation required. On this basis, one would predict the
Type IV matrix to be easier to solve than the Type HI matrix, which in fact is the case.
A third resul suggested by Table 1 is that the amount of computation required to achieve a given
MPPE increases approximately as N 2 . There are approximately N 2 elementary computations (multiply
and add being counted as an "elementary computation") per iteration, so that the amount of computa
tion required in n iterations is proportional to nN 2 . As is evident in Table 1, for a given combination of
matrix type, Po and a, n increases slowly and at a decreasing rate, or does not increase at all, as /Vis
increased. The worst case from the point of view of the relative increase in computation with N is the
Type IV matrix with a = 0.95 and po = 0.9. However, even in this instance only nine more iterations
are required to achieve 0.5percent error as N is increased by a factor of 10 from 50 to 500. Note, how
ever, that in Table 1 with N — 10, there are a good number of instances where more than 10 iterations
and hence more than N 3 elementary computations are required. Thus the relative advantage of the
method as against conventional method requiring approximately /V 3 elementary computations, would
appear to hold only for large values of N.
Finally it is to be noted that the suggested method seems to be insensitive in an appropriate way
to near zero entries in the transition matrices. This is perhaps a major reason for the rapid convergence.
Type I matrices for which 0.005percent accuracy was achieved in four iterations, particularly illustrate
this result as it is obvious that for N = 500 or N = 1 ,000 the entries are extremely small.
COMPUTING RETURN IN MARKOV CHAIN 241
TABLE I. Iterations Required to Achieve 0.5 and 0.005 percent Maximum Possible Error
Matrix Type
I
II
III
IV
\ N
a \
10
50
100
500
1000
Prob
Po
\ N
a \
10
50
100
500
1000
\ N
\
10
50
100
500
1000
\ N
a \
50
100
500
1000
0.995
7
4
5
4
5
3
4
3
4
3
0.9
0.995
4
3
4
3
4
3
4
3
0.995
83
46
90
51
91
51
0.995
58
30
61
32
49
32
0.99
7
4
5
4
5
3
4
3
4
3
0.99
4
3
4
3
4
3
4
3
0.99
80
45
87
49
87
49
0.99
43
19
54
24
0.95
7
4
5
4
5
3
4
3
4
3
0.95
4
3
4
3
4
3
4
3
3
3
0.95
65
37
66
38
69
39
0.95
36
16
43
20
37
25
0.90
6
4
5
4
5
3
4
3
4
3
0.90
4
3
4
3
4
3
3
3
3
3
0.90
50
28
52
30
52
29
0.90
25
12
37
17
29
20
0.5
0.995
5
3
5
3
5
3
6
3
0.995
15
9
16
9
16
9
0.995
15
9
19
10
0.99
5
3
5
3
5
3
5
3
0.99
15
9
16
9
16
9
0.99
11
6
12
7
0.95
5
3
5
3
5
3
5
3
4
3
0.95
15
8
15
9
16
9
0.95
11
6
12
7
12
7
9
0.90
5
3
5
3
5
3
4
3
4
3
0.90
15
8
15
9
14
9
0.90
10
6
11
6
12
7
8
0.1
0.995
7
4
6
3
6
3
5
3
0.995
7
4
6
4
6
4
0.995
6
4
6
3
0.99
7
4
6
3
6
3
5
3
0.99
7
4
6
4
6
4
0.99
6
4
5
3
0.95
7
4
5
3
5
3
5
3
4
3
0.95
7
4
6
4
6
4
0.95
6
4
5
3
5
3
4
0.90
7
4
5
3
5
3
4
3
3
3
0.90
7
4
6
4
6
4
0.90
5
3
5
3
5
3
4
4
4
3
REFERENCES
[1] Kunz, Kaiser S., Numerical Analysis (McGrawHill, New York, 1957).
[2] MacQueen, James B., "A Modified Dynamic Programming Method for Markovian Decision Problems," J. Math. Anal, and
Appl. 14, 3843 (1966).
[3] Scarf, Herbert E., "A Survey of Analytical Techniques in Inventory Theory." In Multistage Inventory Models and Techniques,
Herbert E. Scarf, Dorothy M. Gilford, and Maynard W. Shelly, editors (Stanford Univ. Press, 1963).
[4] Wagner, Harvey M., Michael O'Hagen, and Bertil Lundt, "An Emperical Study of Exactly and Approximately Optimal In
ventory Policies," Management Science, 11, 690723 (1965).
MORE ADO ABOUT ECONOMIC ORDER QUANTITIES (EOQ)
Victor J. Presutti, Jr.
and
Richard C. Trepp
Operations Analysis Office
Headquarters, Air Force Logistics Command
ABSTRACT
This paper is concerned with the determination of explicit expressions for economic
order quantities and reorder levels, such that the cost of ordering and holding inventory is
minimized for specific backorder constraints. Holding costs are applied either to inventory
position or onhand inventory, and the backorder constraint is considered in terms of the total
number of backorders per year or the average number of backorders at any point in time.
Through the substitution of a new probability density function in place of the normal p.d.f.,
explicit expressions are determined for the economic order quantities and the reorder
points. The resulting economic order quantities are independent of all backorder constraints.
It is also concluded that under certain conditions, the minimization of ordering costs and
inventory holding costs (applied to inventory position), subject to a backorder constraint,
is equivalent in terms of reorder levels to minimization of the safety level dollar investment
subject to the same backorder constraint.
INTRODUCTION
We are concerned with a singleechelon, multiitem, continuous review inventory system in which
the process generating demands does not change with time. When an item's reorder level is reached,
an order is placed. After the elapse of a leadtime, which is assumed constant, the entire order is re
ceived. For each item, it is desired to determine the constant reorder quantity and reorder point which
will minimize (summing all items in the inventory) the holding and procurement costs subject to an
overall backorder constraint. The holding cost may be a function of either onhand inventory or inven
tory position. We will constrain total units backordered or average number of units in a backorder
position, which, as is well known, is also the unityears backordered per year. In either case, the
backorders may be weighted by item essentiality.
MODEL DERIVATION
Definitions
The following definitions apply throughout the discussion:
Q: order quantity (units)
r: reorder level (units)
D: annual demands (units)
k: safety factor
cr: standard deviation of unit demand during a leadtime
a: holding cost factor
A: ordering cost (dollars per order)
c: item cost (dollars per unit)
fju: mean leadtime demands
Z: item essentialtiy (relative military worth)
243
244
V. J. PRESUTTI, Jr. AND R. C. TREPP
x: units demanded during leadtime
T)\ constraint — total essentialityweighted units backordered
/3: constraint — expected number of essentialityweighted units in a backordered position
at any point in time
t: leadtime.
We also require:
i) r = fx + kcr
ii)/(x),the probability density function (p.d.f.) describing the distribution of units demanded
during the leadtime:
(1)
2^ ex P (
X — /JL
for all x.
It is easily shown that this implies:
(2)
and
(3)
Pr[x > fx+ ko]= 0.5 exp ( V2k) for k 2*
Pr[x > fJL + k<r] = 1  0.5 exp ( V2A) for k < 0.
The form of f(x) was inspired by the work of L. L. Parker [3] who essentially assumed demands
were normally distributed, and then approximated the complementary cumulative function by a function
of the form ae~ bx . We originally chose a and b such that the mean absolute deviation from the normal
distribution over the range =s k =s 3 was minimized (a = 0.57, 6=1.35). We now consider, however,
that a = 0.5 and b— v2 is more appropriate.* In addition to defining a p.d.f. that is symmetric about
the mean (/u,), these values also cause Pr
hi < ?—£ < k 2
to be independent of fi and u (a desirable
attribute also shared by the normal p.d.f.). Further, with these values, the expected backorders at a
random point in time for k = is the same as if demands were normally distributed, for large Q/cr. This
similarity was considered desirable since we are primarily interested in average backorders (rather than
total backorders) and we will often have no safety level.
The value of P r \k =£ =£ k + 0.2 for k = to 4.0 is given in Table 1 for both the normal p.d.f.
and f{x).
We will now operate as if
A.)^(vi(^)r
*We are indebted to Mr. Alan Kaplan of the AMC Inventory Research Office who suggested that we reconsider our values
for a and b.
**The absolute value of — ■ — is omitted.
EOQ
245
Table 1
Comparison o/f(x) to the Normal p.d.f.
K
A"
B"
C n (xY
C f (xY
0.0793
0.1232
0.0793
0.1232
0.2
0.0761
0.0928
0.1554
0.2160
0.4
0.0703
0.0700
0.2257
0.2860
0.6
0.0624
0.0527
0.2881
0.3387
0.8
0.0532
0.0397
0.3413
0.3784
1.0
0.0436
0.0299
0.3849
0.4084
1.2
0.0343
0.0226
0.4192
0.4310
1.4
0.0260
0.0170
0.4452
0.4480
1.6
0.0189
0.0128
0.4641
0.4610
1.8
0.0132
0.0097
0.4773
0.4704
2.0
0.0088
0.0073
0.4861
0.4777
2.2
0.0057
0.0055
0.4918
0.4832
2.4
0.0035
0.0041
0.4953
0.4874
2.6
0.0021
0.0031
0.4974
0.4905
2.8
0.0013
0.0024
0.4987
0.4928
3.0
0.0006
0.0018
0.4993
0.4946
3.2
0.0004
0.0013
0.4997
0.4959
3.4
0.0001
0.0010
0.4998
0.4969
3.6
0.0001
0.0008
0.4999
0.4977
3.8
0.0001
0.0006
0.5000
0.4982
4.0
0.0000
0.0004
0.5000
0.4987
— =S/t + 0.2 \fom(x)=— i=exp
•V2^
\ 2<r 2 /
■B = Pr\k* ^k+0.2\forf(x)=— c .exp(V2
cr J 2cr \ cr
x — fj.
xV
}
X/A
C(x)=Pr =£ s ArH0.2l for n(x) and /(*), respectively.
246 V. J. PRESUTTI, Jr. AND R. C. TREPP
Notice that this is strongly biased against negative k. In fact, by operating as if the probability dis
tribution could be described by Eq. (4), we are, in effect, saying that limit Pr[x > p. + kcr] = °°, a
k— *— °°
somewhat preposterous allegation. Nonetheless, this bias is considered desirable since in actual
practice, we restrict k 3 s 0.*
Expected Backorders
In order to determine the expected number of units backordered per year (Bx) and the average
number of units in a backorder status at a random point in time (Bt) , we need g(y) —the p.d.f. for
the number of backorders at a random point in time. The following technique for obtaining g(y) was
developed by Hadley and Whitin [1].
g(y)dy— probability that the number of units backordered is between y and y + dy, and
£(Vl(«±5=£))*
= the probability that the number of units backordered at time, t, is between y, y+ dy given that
the inventory position (onhand plus duein minus dueout) is between p, p + dp at time t — t. But the
probability that the inventory position is between p, p + dp is dpIQ.
Therefore,
^r^e, P (V2(^)), r *
and
and finally
(7) g(y) = ^ exp ( V2 ( Z ^)) d " exp( V2QI<r)).
Consequently, the probability of a stockout (P ou t) is
(8)  o " g(y)dy = (^)  (1  exp ( V2QI*)) exp ( V2k).
Therefore,
(9) B x = DP 0Ut = (~^j Q (1  exp ( V2QI<t) ) exp ( V2k)
(seeRef. [1], page 178)
and
(10) B T = f~ yg{y)dy = ^^j (1exp ( V2QI<r)) exp ( V2A).
Other Inputs
We must now obtain the average number of procurements per year, the expected inventory
position (E[IP]), and the expected onhand assets (E[OH]). Clearly, the average number of procure
ments is D/Q.
*The nonnegativity requirement for A: is a result of the desire to make this model acceptable to practicing logisticians.
EOQ
247
It is also well known that
and
By definition,
Therefore,
£[Duein] = \x (see Ref. [1], page 187)
E[IP]=r+QI2 = ix + ka + QI2.
E[IP] =E[OH] +£[Duein] E[B/0].
E[OH] =E[IP] £[Duein] +E[B/0]
= fi + ko + QI2ix + B T
= k<r + QI2 + B T .
FOUR INVENTORY MODELS
We are now ready to proceed to the four models which will enable us to determine Q and r for each
item in such a way as to minimize the holding cost (applied either to E[OH] or E[IP]) plus the pro
curement cost, subject to an overall backorder constraint. This constraint may be either the total units
backordered per year (i.e., the stockout penalty is independent of the duration of the stockout) or the
average number of units in a backorder position (i.e., the stockout penalty depends on the duration of
the stockout). In either case, the backorders may also be weighted by item essentiality.
MODEL I: Stockout penalty independent of duration of stockout; holding cost applied to onhand
inventory. We wish to minimize
(LI)
subject to
(1.2)
A AiD t 
ZnT + Z aiCi
i=l v ' t=l
kcri+y+^^(lexp(V2^)
A 0.5 ZiDjOj
(lexp(V2^) exp ( V2ki) *S t,
£V2 Q t
where Z, represents item essentiality. The Method of Lagrange yields
QidiCi
(1.3)
(1.4)
*' = "V? ln
V^))(^pXZ (J D,
2 (afim =V27,,and
VJ
kZiDi
(1.5) o = ^ + ^i
Qf 2
aidCTi
V2" 1exp [V2
£))
1exp [V2
aiJJ V2
Qi
<Ti
exp
(^5)
For reasons which will soon become evident, we will not solve (1.5) for Q t until later.
MODEL II: Stockout penalty independent of duration of stockout; holding cost applied to in
ventory position.
"AM" ( Q t
Z ~q~ + Zj aiCi[fJ<i + ki(Ti + Y
248 V. J. PRESUTTI, Jr. AND R. C. TREPP
We wish to minimize
(III)
subject to
" 0.5 ZiDm I ( _ Qi\ ,
The Method of Lagrange yields
1
(113)
"v^ ln
QidjCj
0.5 (A) Z,Z) ( flexp ( V2"
<?.■
(H.4)
A = — — — V cnatCi, and
V2 Ti £i
i"
(115) = ^i + ^
aiCicn
Vf(lexp(V2a))
1 — exp
V2&
MODEL III: Stockout penalty depends on duration of stockout; holding cost applied to onhand
inventory.
We wish to minimize
(IH.l)
subject to
(III.2)
)^+;aiC i [fe(ri+^+^g(lexp(V2^))exp(V2A i )]
i^^(lex P (V2^))ex P (V2fe)^ i 8.
The Method of Lagrange yields
(III.3)
kt — ;= In
V2
V2(W,
0.5(7, ( 1 exp ( V2 & ) ) (aidkZi)
(III.4)
1,(5^=^. and
(III.5) °=^ +S f
aiCicn
V2(l
exp
V2
CTi
1 — exp
V2&
<?,
^expfV^
oj
CXi
MODEL IV: Stockout penalty depends on duration of stockout; holding cost applied to inventory
position.
EOQ
249
We wish to minimize
(IV.l)
subject to
(IV.2)
y jt + i aiCi ( /*■• + kiai + 9
i=l v* i = 1 \ ^
The Method of Lagrange yields
(IV.3)
ki= ~vz ln
V2 QidiCi
0.5 ( \)ZiO, 1 exp  V2
(J,
(IV.4)
(IV.5)
" cr ( aiCi
: > — ;^t:, and
<?? 2 "
atCjCTi
V? (lexp  V2"
a
<Ti
lexp lv^
5i
^exp(
(3i
cr,
Notice that in all four models Qi is the same (Eqs. (1.5), (II. 5), (HI. 5), and (IV.5)). The quantity
Qi is completely independent of the backorder constraint regardless of whether or not the constraint
is timeweighted. We also note that we cannot solve explicitly for Q. Conceptually, there is no problem.
The exact solution for Q can be obtained through the use of an iterative routine (assuming access
to a computer). An alternate and perhaps a more suitable approach is to recognize that in most cases
(! p(Vjfi))
]
and can therefore be ignored. Consequently,
AiDi
(13)
and
(14)
atCi ajCiCTi
~ ~ Qf + 2 V2C,
Qt V2 + V a iCi + lv2T
Notice that as <t, — > 0, Eq. (14) approaches the well known Wilson lot size formula.
DISCUSSION MODEL APPLICATION
Holding Cost Considerations
We wish now to discuss the difference between applying holding cost to the onhand inventory
(Models I and III) and applying it to the inventory position (Models II and IV). In spite of tradition, we
feel it is more logical to apply the holding cost to the inventory position (see Ref. [4] for an argument
supporting this point of view). From a practical point of view, the difference in results is small; however,
the inventory position models have two very desirable attributes. First, there is an explicit solution
*It is easily shown that this value is an upper bound for the exact solution.
250
V. J. PRESUTTI, Jr. AND R. C. TREPP
for X, the Lagrangian, in these models — which is not the case when the holding cost is applied to the
onhand inventory. Second, in the inventory position models the safety level expressions are equivalent
to those that would have been obtained had we somehow arrived at the buy quantities ((Vs) and then
independently minimized the cost of holding safety stock subject to the same backorder constraints
used in Models I through IV. That is, under the inventory position models, the inventory manager
who minimizes backorders subject to a safety level cost constraint would arrive at A:,'s equal to the A;,'s
that he would have calculated had he minimized the total costs of ordering and holding inventory subject
to a backorder constraint. This is especially important when the optimal (Vs are not used: for example,
when the policy is to purchase 3, 6, 9, or 12 months' worth of stock, or even when Qi is determined by
the Wilson lot size formula. If the same holding cost factor is applied to all items (i.e., at = a for all t),
the results obtained in models II and IV are equivalent to those determined when minimizing the safety
level investment subject to a backorder constraint.
Negative k's
Although all four models are strongly biased against negative &'s, as long as there is a heterogeneous
inventory with some high cost, low demand items, there will be some items for which the models will
specify negative &'s. It could be argued that if even with a strong bias negative A:'s are obtained they
should be accepted. We have no philosophical quarrel with this argument; however, in this paper, we
insist that there be no negative safety levels. In order to accomplish this, we will set all negative A's
equal to zero. Consequently, we will have more protection than the models call for and thus have fewer
backorders, but a higher cost. The problem is to determine the X that will allow, taking into account the
nonnegativity of k, an acceptable number of backorders at a reasonable cost. Since we believe that
Model IV is most appropriate for the Air Force's nonreparable items, we will use it to describe the
methodology for determining A.
Employing a simulation using actual data for about 5,000 Air Force items (for which seven years
of history is available) we will select a range of values for [—A.]. Each value of [—A] will yield 5,000
A;'s — one for each item. The nonnegative A:'s will be accepted and the negative &'s will be set equal to
zero. We will then, using these A;'s which are now all nonnegative, obtain the average number of
essentialityweighted units in a backorder status divided by the number of essentialityweighted units
demanded annually (call this ratio /?«), and the cost of safety stock divided by the dollar value of annual
demands (call this ratio R t ). Subsequently, Figures 1 and 2 will be plotted as shown below:
3t
Given Figures 1 and 2, the inventory managers can select either R B or R$ (hopefully they will make
the selection by considering both the cost and the effectiveness aspects) and then choose from the
appropriate graph the value of [— X]. The ability to choose on the basis of fund availability is extremely
important since, in practice, the inventory manager is considerably more likely to be subject to a
system constraint of dollars available for safety level investment (out of pocket cost) than to a back
eoq 251
order constraint. Hadley and Whitin emphasize this point in their* discussion of budget constraints.
They conclude:
"None of the inventory models discussed in operations research literature appears to provide
an adequate treatment of budgetary constraint. In the first place, budgets are usually set in
terms of actual dollars to be spent, whereas the costs included in inventory control models
are often not the outofpocket cost type."*
An example of how the inventory manager (IM) might use these graphs (possibly to reduce "actual
dollars to be spent") is as follows:
Suppose the IM submits a budget request, for operating and safety stock, which allows 20 percent
of the dollar value of annual demands for safety stock. If he receives the monies requested he can
obtain from Figure 2, the [— A.] that will yield an Rt of 20 percent. If he gets M dollars less than
requested, and that M dollars represents 5 percent of annual dollar demands, he can set R%
equal to 15 percent (20% — 5%) and obtain the new value for [— \].
While this paper doesn't concern itself (other than the portion immediately preceding) with
constrained optimization, we have applied the demand distribution described in this paper to models
that do constrain the buy quantities (i.e., the (?'s), the onhand inventory, and the inventory position
as well as the safety stock. None of these methods, however, is satisfactory when it comes to reducing
actual dollars to be spent by very large amounts below computed requirements. All of these constraining
models reduce dollars spent by reducing the reorder level (i.e.. by reducing k) and/or by reducing the
order quantities in some manner. We cannot, however, reduce the reorder level greatly without paying
a very large price in terms of backorders. If most of our expenditures were for items with multiyear
buy quantities, then we could temporarily reduce actual dollars to be spent by reducing these buys
to a year or less. Unfortunately, by far the greatest majority of dollars are spent on items for which the
buy quantity is on the order of three months worth of stock — which theoretically eliminates the pos
sibility of obtaining a large reduction in actual dollars to be spent by reducing the buy quantity.
Backorder Costs Known
The reader will note that in the event the stockout cost (s,) is known, the problem is simplified.
The (?,'s remain unchanged and the A,'s of each of the models have a 5; in place of [— A] Z, (see Eqs.
(1.3), (II.3), (III.3), and (IV.3)). There is no [ A] .
Final Note
The Air Force is planning to test all four models, nevertheless, at the present time, it favors (un
officially) Model IV. This model has also been brought to the attention of the Department of Defense
by the Air Force for consideration during a series of recent joint services meetings designed to further
standardize inventory policies and objectives among all DoD components.
REFERENCES
[1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, New Jersey, 1963).
[2] Hadley, G. and T. M. Whitin, "A Review of Alternative Approaches to Inventory Theory," The RAND Corporation, Santa
Monica. California. RM4185PR (Sept. 1964).
(3 Parker, L. L.. "Economical Reorder Quantities and Reorder Points with Uncertain Demand," Nav. Res. Log. Quart.. //,
351358 (Dec. 1964).
[4] Schrady, D. A., "Continuous Review Inventory Practice in the Navy Application," Report #NPS55o9031A, United States
Naval Postgraduate School, Monterey, California (Mar. 1969).
* Reference [2, page 152.
THE RELATIONSHIP BETWEEN DECISION VARIABLES AND PENALTY
COST PARAMETER IN (Q, R) INVENTORY MODELS
Alan J. Kaplan
AMC Inventory Research Office
Army Logistics Management Center
Fort Lee, Virginia
ABSTRACT
This paper is concerned with the optimum decision variables found using order quantity,
reorder point (Q, R) inventory models. It examines whether the optimum variables
(Q* and R*) are necessarily monotonic functions of the backorder cost parameter (or
equivalently of the performance objective). For a general class of models it is proved that
R* must increase as the performance objective is raised, and an inequality condition is
derived which governs how Q* will change. Probability distributions of lead time demand
are cited or found for which Q* increases, Q* decreases, and Q* is independent of increases
in performance objectives or backorder cost parameter.
INTRODUCTION
This paper is concerned with the optimum decision variables found by using order quantity,
reorder point (Q, R) inventory models. It examines whether they are necessarily monotonic functions
of the penalty (backorder) cost parameter (or equivalently of the performance objective). It attempts
to make general statements about the direction of such monotonicity. Such information is of particular
interest in contexts where the true cost of a backorder is difficult to estimate — and the performance
objective is determined in large part by its impact on inventory levels.
It is easy to be misled by intuition. Our own empirical experience was with a (Q,R) model in which
backorder cost is proportional to the length of time a backorder is outstanding, and lead time demand
is normally distributed. For this case, the optimum reorder point (R*) always increased as the penalty
cost parameter (P) increased, while the optimum order quantity (Q*) decreased. This is illustrated
in Herron's [2] dimensionless graphs for the normal distribution. The conclusion can be proved analyti
cally for deterministic demand and penalty proportional to time on backorder by solving directly for Q *;
e.g., see Hadley and Whitin's [1] solution.
On the other hand, Victor Presutti [3], found Q* to be independent of the penalty cost parameter
when he used a pseudoexponential distribution for lead time demand. t That finding prompted this
paper.
We initially assume penalty cost is proportional to the length of time a backorder is outstanding
and begin by reviewing the relevant aspects of the appropriate (Q, R) model. Next we obtain some
useful theoretical conclusions, including a proof that R * must always increase when P does and an
inequality condition which governs how Q* will behave. We apply these results, and then in the last
section extend them to other {Q, R) models.
t/(x) = A.ie x " forx^O, and/(jc) denned elsewhere so that / f(x)dx=\. Presutti's work is based on Lawrence Parker's
approximation to the normal complementary cumulative.
253
254 A  J KAPLAN
The theoretical conclusions require the use of an approximation for average backorders (Eq. (2)).
This approximation is commonly used in the actual computation of /?* and Q* [1], [2J, and does not
negate the value of the results.
THE MODEL
The model is the basic steady state, continuous review (Q, R) model applicable when penalty
cost is proportional to the length of time a backorder. is outstanding. For the model to give strictly
optimal answers in the general case, procurement lead times must be constant and demands must
be for one unit at a time. However, the model is being used profitably even when these conditions do
not hold. For ease of exposition we will assume that the lead time demand distribution is continuous.
One good reference, particularly for providing the basis for Eqs. (1) and (2), is Hadley and Whitin [1].
Let
TC = average total cost per year;
H= holding cost per year per unit of stock on hand;
P= penalty cost per year per unit of stock on backorder;
A = administrative and/or setup cost;
A — average yearly demand;
A/. = expected demand during procurement lead time;
B — average amount on backorder at a random point in time;
G(x) — probability that lead time demand equals or exceeds x;
L = lead time.
Then
TC(R,Q)=A(klQ)+H (R + QI2\,. + B) +P • B = A • {KIQ)+H (R + QI2k,.)+ P'B
(1)
P' = P + H
B = ±r (xR)[G(x)G(x + Q)]dx.
v Jk
For computational purposes, the following approximation is often used for B:
(2 ) B = jjf~(xR)G(x)dx.
The approximation is valid to the extent that the probability of demand exceeding R + Q is small.
We will utilize this expression for B.
THEORETICAL RESULTS
Let
Qi, R, be the optimal decision variables for P'=/ > t ;
Qi, Ri be the optimal decision variables for P' =P 2 ;
Pz>Pr,
Bt be expected average backorders for decision variables Qi, Rr,
F(Q,R)
(Q,R) INVENTORY MODELS 255
dB/dQ
SB/dR'
The following two theorems will be proved:
THEOREM 1: R* (the optimal R) must increase as P increases.
THEOREM 2: The sign ofdF(Q, R)/dR, if it is consistent, determines whether ()* will also increase.
The key to these results is the corollary to Lemma 1 that the algebraic sign of Q 2 — Q\ is determined
by the algebraic sign of F(Q 2 , R 2 ) — F(Q U R t ). Lemma 2 and its corollary demonstrate the intuitively
"obvious": that expected backorders (B) must decrease in the optimum as P increases and, therefore,
that it is not possible for both Q* and R* to decrease. Lemma 3 shows that dF{Q, R)ldQ is always
negative. Lemma 2 is used in proving Theorem 1, while Lemma 3 is used in proving Theorem 2.
LEMMA 1: The algebraic sign of F(Q 2 , Rz) — F{Q\, R\) is the same as the algebraic sign of
Qi~Q\ (A zero value for one implies a zero value for the other.)
PROOF: Taking the partial derivatives of the total cost equation, Eq. (1), with respect to Q* and
/?*, and setting the derivatives to zero gives two conditions:
(3a) (dBldQ*)P'^AI(Q*y + H/2 = 0;
(3b) (dB/dR*)P' + H=0.
By algebraic manipulation (transposing and dividing one equation by the other), this gives
F(Q*,R*) = [A AI(Q*) 2 HI2] +H.
Hence
F(Q 2 , R 2 ) F(Q U /?,) = [XAKQ.rXAKQ.y] //.
If Q 2 exceeds Qi(i.e., Q2 — Q1 is positive), then both the numerator and denominator on the right
hand side of Eq. (4) are negative; F(Q 2 , R 2 ) — F(Qi, Ri) is then positive. By similar reasoning the
rest of Lemma 1 is shown to hold.
COROLLARY: The algebraic sign of Qt — Q\ is the same as the algebraic sign of F(Q 2 , R 2 )
— F(Q\, Ri). A zero value for one implies a zero value for the other.
LEMMA 2: The expected value of backorders (in the optimum) must decrease as P increases.
PROOF: To simplify the explanation let Oi stand for operating costs corresponding to decision
variables (Q u Ri) so that Total Cost = O t + B t ■ P' .
We use proof by contradiction. We initially assume B 2 > B t . By the optimality of Q u Ri for/ > '=P 1 ,
Using P, > P, and B>>Bi, we find
O. + BrP^BriPiPt) > O. + fi.P. + fi.^PiP,).
Hence
2 + B 2 P 2 >0 1 + Bi P 2 .
But this last inequality contradicts the optimality of Q 2 , R 2 when P' =P 2 .
COROLLARY: Both Q * and R * cannot decrease as P increases (proved by Lemma 2 and inspection
of Eq. (2)).
LEMMA 3: dF{Q, R)/dQ is always negative.
PROOF: Using Eq. (2), we get by calculus,
256 A. J. KAPLAN
(5a) % = ~7p\l ('WW**
dB 1 f v ,
= ~Q in G(x)dX 
and
(5b)
Hence,
dR
(6) F(Q, R)
~j^(xR)G(x)dx ±j^ (xR)G(x)dx g
— I™ G{x)dx !" G(x)dx [ X G{x)dx
(JjR JR JR
The denominator does not depend on Q. Clearly, the numerator decreases as Q increases.
THEOREM 1: R* increases as P increases.
PROOF: Using Eq. 6, we have
(7a) F(Q 2 ,R 2 )=B 2 +\G(x)dx
J «2
and
(7b) F(Q u R 1 )=Bi^r('G(x)dx.
JRi
Suppose R 2 < Ri. Then the denominator of (7a) will be larger than the denominator of (7b). By
Lemma 2, B 2 <B\, so this means F{Q 2 , R 2 ) <F(Qi, Ri). By the corollary to Lemma 1 this implies
Qz<Q\. By the corollary to Lemma 2, however, it is not possible for Q 2 < Q\ if R 2 < Ri since P 2 > Pi.
Proof by contradiction.
THEOREM 2A: If dF(Q, R)/dR is consistently negative, then Q* must decrease as P increases.
Since R 2 >Ru and dF(Q, R)/dR is given as negative, F(Q U R 2 ) <F(Q l , R^. Now if Q 2 ><?,,
then by Lemma 3 this means F(Q 2 ,R 2 ) < F(Qi, R 2 ). By combining the two inequalities, we get
F(Q 2 , R 2 ) < F(Qi, Ri). However, by the Corollary to Lemma 1, this implies Q\>Q 2 . Proof by
contradiction.
THEOREM 2B: If dF(Q, R)/dR is consistently positive, then Q* increases with P.
The proof is analagous to that of Theorem 2 A. Since R 2 > Ri, and dF(Q, R)ldR is given as positive,
F(Q U R 2 ) >F(Q U /?,). Now if Q 2 <Q U then by Lemma 3 this means F{Q 2 ,R 2 ) >F(Q U R 2 ). By
combining the two inequalities we get F(Q 2 , R 2 ) > F(Q t , R\). But by the Corollary to Lemma 1 this
implies @i < (?2. Proof by contradiction.
THEOREM 2C: If F{Q, R) is independent of R, then Q* is independent of P. For this case
F(Q U R 2 )=F(Qi, Rt). If <? 2 < <?i then by Lemma 3, F((? 2 , /? 2 ) > F{Q X ,/?,), which by the Corollary
to Lemma 1 implies Q 2 >Q\. Similar reasoning shows that the assumption ()i < Q 2 also leads to a
contradiction.
Applications
There is a strong temptation to reason intuitively that Q* must increase as penalty cost increases:
higher reorder quantities reduce the probability of stock out; this advantage becomes more significant
as the penalty cost parameter is increased, while the impact of Q on ordering and holding costs is
unchanged. Why is this reasoning wrong?
The flaw is that when penalty cost increases, it becomes worthwhile to increase /?*, so that chang
ing Q* has less impact on backorders than before, and therefore may have less impact on penalty
(Q,R) INVENTORY MODELS 257
costs (the product of backorders and the penalty cost parameter). The decisive factor relates to the
relative effectiveness in reducing backorders of augmenting Q* vs augmenting R*. As R increases,
does the relative effectiveness of augmenting Q increase? Only then will Q* increase as P increases.
It is difficult to apply Theorem 2 analytically to most probability distributions, i.e., to determine
analytically the sign (+, — , or mixed) of dF(Q, R)dR. No attempt was made to carry on an extensive
computer aided mapping of sign vs probability distribution
We did attempt to find a case when Q* would increase as penalty cost increases. Recall that for
the normal distribution Q* decreases and that for exponential type distributions Q* is unchanged.
Let the Weibull distribution be defined:
f(x) = aXx" 1 exp ( Xx a ); x^0,Xa^0.
Then for a= V2, we have analytically shown that dF(Q, R)/dR is positive for R 3= 0, i.e., this is a
case where Q* increases with R. Empirical work indicates that this conclusion about the Weibull
distribution also holds true when the precise expression for average backorders (Eq. (1)) is used.t
For a=l the Weibull reduces to the exponential. It is simple to show dF(Q, R)/dR — for the
exponential regardless of which expression for B is used.
As a matter of interest, the Weibull with a= V2 was originally examined because of an inference
made from an analogy between determination of the sign of dF(Q, R)/dR, and the sign of the deriva
tive of the failure rate function of reliability theory, the latter being a standard problem.
EXTENSIONS
The general linear form for penalty cost is P + P • t where t is the time on backorder. This trans
lates into total expected penalty cost =P ■ B + P ■ B or, in words (cost per unit backordered per year
on backorder) X (average amount on backorder) + (fixed cost per backorder occasion) X (number of
demands backordered each year). That part of the total cost equation, Eq. (1), dependent on B or B
becomes P' B + P ■ B where P' = P + H.
The theory derived in this paper still holds provided that as the penalty cost parameters are
increased, P remains equal to k • P' for constant k (or P — 0, the special case already considered).
If the cost item H • B is left out of the total cost equation, Eq. (1), the condition is simply that P — k • P,
which is more natural. H • B is normally very small, and in fact does not appear as a cost if holding
cost depends on inventory position (including net duein) rather than on hand.
An approximate expression must be used for B as well as B. 1 1
Let
fjl'W
B T =kB + B.
For F T {Q, R), defined in terms of dBrldQ and BB T lBR, we have
Br
Ft(Q,R) =
(*G(x)dx + kkG(R)
tMiss ChungMei Ho did much of the analytical investigation on the application of Theorem 2, and Mr. Robert Deemer
did the empirical work for the Weibull.
ttHadley and Whitin, chapter 4 has a derivation.
258 A. J. KAPLAN
By inspection, all three lemmas proved hold with the substitutions of Bt for B, P' • Bt for P'B, and
F T (Q, R) for F(Q, R). Therefore, the two theorems hold with these substitutions.
Using the precise expression for B causes difficulties even when P is 0. Then
BT [y(R + Q)]g(y)dy , C( .
F{Q,R)=—±** , where g ( y ) =&&.
[G(x)G(x + Q)]dx . y
It is no longer apparent that Lemma 3 must always hold, or that in the proof of Theorem 1 that
the numeraton F(Q 2 , Rz) < numerator a£F{Qu Ri) if Qi > Q\
REFERENCES
[1] Hadley, G. and T. M. Whitin, Analysis of Inventory Systems (Prentice Hall, Englewood Cliffs, N.J., 1963).
[2] Herron, David P., "Use of Dimensionless Ratios to Determine MinimumCost Inventory Quantities," Nav. Res. Log. Quart.
13, 167176(1966).
[3] Presutti, V. J. and R. C. Trepp "More Ado about Economic Order Quantities (EOQ)", Nav. Res. Log. Quart. 17, 242250
(1970).
NEWS AND MEMORANDA
PROFESSIONAL DESIGNATION IN LOGISTICS MANAGEMENT
The School of Systems and Logistics, Air Force Institute of Technology, in cooperation with
the Society of Logistics Engineers (SOLE), has announced the beginning of a program of professional
designation for qualified military and Civil Service logistics managers of the Department of Defense.
Those who meet the criteria may be awarded an appropriate certificate of "Professional Designation
in Logistics Management." Certification is not restricted only to members of the Society.
The purpose of the program is to recognize the efforts of the individual to increase his educa
tional background in the broad world of logistics. It recognizes that true logistics management capa
bility must encompass more than expertise in one of the several disciplines which comprise the logistics
spectrum. The idea is that logistics is inherently an integrating process of these several disciplines
and that the effective logistics manager will be one who has knowledge in, and understands the view
point of, the various subfunctions which the separate disciplines represent. The professional stature
of the functional specialists is recognized, but this program is aimed at the need lor a portion of these
managers to become logistics generalists to perform the essential function of integrating the efforts
of many functional specialists.
Qualification for certification rests upon the types and forms of education and experience which
contribute to the overall concept of integrative management in logistics. The courses of study outlined
in the qualification listing were tailored with this concept in mind. The objective is to assure that those
who qualify have progressed from one of the functional specialties to a broader realization of logistics
management needs through exposure to education in related and supporting disciplines and functional
areas. Thus, those who qualify for certification will have demonstrated a high degree of initiative in
acquiring a broad educational and experience background. They will be considered to have demon
strated their competency for advancement to higher levels of responsibility in logistics management
where their broadened point of view would aid decisionmaking of more costeffective nature.
The certification program provides for two levels of certificate award. The first is "Professional
Designation in Logistics Management" and the second is "Professional Designation in Logistics
Management — Advanced." The qualification criteria includes courses in five categories as indicated:
\ Certificate Advanced
Category A Logistics Integration 2 2
Category Bl Materiel Acquisition 1 2
Category B2 Materiel Distribution 2 2
Category B3 Materiel Maintenance 1 2
Category C Management Techniques 2 2
Total Required 8 10
Some allowance is made in the program for waiving certain educational course requirements on
the basis of field experience. Such waiver, for a maximum of two courses, must be substantiated by
the applicant. The determination of whether or not to accept the experience equivalency is left to the
discretion of the certification committee of the School of Systems and Logistics. Applicants are urged
259
260 NEWS AND MEMORANDA
to explain their request for equivalency and waiver by an adequate explanation of the experience
(where, for how long, in what job, at what level, and so on).
It is also likely that many will have taken higher level educational courses than those listed in
the certificate categories. These may be accepted in lieu of the described courses if the applicant
includes proof of completion along with a brief outline of material covered.
Certificates will be approved by the School of Systems and Logistics and the SOLE National
Education Committee through a certification board process. Applicants, who need not be members
of SOLE but who must be members of DOD, are encouraged to submit their requests for certification
to the School of Systems and Logistics (AFITSL), WPAFB, Ohio 45433. Each request must be
accompanied by the necessary backup data.
CHANGE IN OFFICE OF NAVAL RESEARCH ADDRESS
Effective with the relocation of Department of the Navy activities from the Main Navy Building
in Washington, the address for the Office of Naval Research, including the Naval Research Logistics
Quarterly, is Arlington, Virginia.
Correspondence for the Naval Research Logistics Quarterly should be addressed as follows:
Managing Editor
Naval Research Logistics Quarterly
Office of Naval Research
Arlington, Virginia 22217
■fr U.S. GOVERNMENT PRINTING OFFICE i 1970 OL386049
INFORMATION FOR CONTRIBUTORS
The NAVAL RESEARCH LOGISTICS QUARTERLY is devoted to the dissemination of
scientific information in logistics and will publish research and expository papers, including those
in certain areas of mathematics, statistics, and economics, relevant to the overall effort to improve
the efficiency and effectiveness of logistics operations.
Manuscripts and other items for publication should be sent to The Managing Editor, NAVAL
RESEARCH LOGISTICS QUARTERLY, Office of Naval Research, Arlington, Va. 22217.
Each manuscript which is considered to be suitable material tor the QUARTERLY is sent to one
or more referees.
Manuscripts submitted for publication should be typewritten, doublespaced, and the author
should retain a copy. Refereeing may be expedited if an extra copy of the manuscript is submitted
with the original.
A short abstract (not over 400 words) should accompany each manuscript. This will appear
at the head of the published paper in the QUARTERLY.
There is no authorization for compensation to authors for papers which have been accepted
for publication. Authors will receive 250 reprints of their published papers.
Readers are invited to submit to the Managing Editor items of general interest in the field
of logistics, for possible publication in the NEWS AND MEMORANDA or NOTES sections
of the QUARTERLY.
NAVAL RESEARCH JUNE 1970
LOGISTICS VOL. 17, NO. 2
QUARTERLY NAVSO P1278
CONTENTS
ARTICLES
Page
Markov Chain Analysis of a Situation Where Cannibalization is the Only Repair Activity by A. J.
Rolre
Some Estimates of Reliability Using Interference Theory by M. Mazumdar
The Optimal BurnIn Testing of Repairable Equipment by J. M. Cozzolino
Goal Programming in Econometrics by W. A. Spivey and H. Tamura
i
Production and Employment Scheduling in Multistage Production Systems by C. Haehling von
Lanzenauer
A Mutual PrimalDual Linear Programming Algorithm by M. Y. Harris
Interception in a Network by R. D. Wollmer
The Fixed Charge Problem by D. I. Steinberg
On Computing the Expected Discounted Return in a Markov Chain by D. F. Hitchcock and J. B.
MacQueen
More Ado About Economic Order Quantities (EOQ) by V. J. Presutti, Jr. and R. C. Trepp
The Relationship Between Decision Variables and Penalty Cost Parameter in (Q,R) Inventory
Models by A J. Kaplan
News and Memoranda
OFFICE OF NAVAL RESEARCH
Arlington, Va. 22217